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ABSTRACT 

Results are presented of numerical simulations of normal isolated late type spiral 
galaxies. Specifically the galaxy NGC 628 is used as a template. The method employs 
a TREESPH code including stellar particles, gas particles, cooling and heating of the 
gas, star formation according to a Jeans criterion, and Supernova feedback. A regu- 
lar spiral disc can be generated as an equilibrium situation of two opposing actions. 
On the one hand cooling and dissipation of the gas, on the other hand gas heating 
by the FUV field of young stars and SN mechanical forcing. The disc exhibits small 
and medium scale spiral structure of which the multiplicity increases as a function of 
radius. The theory of swing amplification can explain, both qualitatively and quanti- 
tatively, the emerging spiral structure. In addition, swing amplification predicts that 
the existence of a grand design m — 2 spiral is only possible if the disc is massive. The 
simulations show that the galaxy is then unstable to bar formation, confirming the 
result of Ostriker & Peebles (1973). The occurrence of this bar instability is further 
investigated. A general criterion is derived for the transition between bar stable and 
unstable, depending on disc mass contribution and on disc thickness. It seems that 
bar stability hardly depends on the presence of gas. A detailed quantitative analysis 
is made of the emerging spiral structure and a comparison is made with observa- 
tions. That demonstrates that the structure of the numerical isolated galaxies is not 
as strong and has a larger multiplicity compared to the structure of some exemplary 
real galaxies. It is argued that the suggestion of Kormendy & Norman (1979) holds, 
that a grand design can only be generated by a central bar or by tidal forces resulting 
from an encounter with another galaxy. 

Key words: galaxies: evolution - galaxies: general - galaxies: fundamental parame- 
ters - galaxies: kinematics and dynamics - galaxies: structure 



1 INTRODUCTION 

It is now well established that substantial amounts of dark 
matter are needed to explain flat rotation curves in the outer 
regions of spiral galaxies. However, the amount of dark mat- 
ter present in the inner, luminous region of a galaxy is not 
well determined. Does a spiral galaxy conform to the maxi- 
mum disc hypothesis (van Albada & Sancisi 1986; Salucci et 
al. 1991; Sellwood & Moore 1999) which states that the lu- 
minous part should be scaled up as much as possible? Or is 
the maximum rotational contribution of the disc to the total 
rotation around 63%? The latter being derived from obser- 
vations of disc stellar velocity dispersions (Bottema 1993, 
1997) or suggested by a statistical analysis of rotation curve 
shapes in relation to the compactness of discs (Courteau & 
Rix 1999). When discs obey the 63% criterion, their mass- 
to-light ratio is half that of the maximum disc case, but the 
disc still dominates in the inner two to three scalelengths. 
A rotation curve analysis generally allows a maximum disc 
rotational contribution anywhere between 50 to 90% (van 



der Kruit 1995). The above mentioned velocity dispersion 
measurements have only been done for a dozen galaxies and 
the analysis depends on the not well determined ratio of 
disc scalelength to scaleheight. As observations are not yet 
conclusive, numerical simulations of spiral galaxies may give 
clues regarding the disc to halo mass ratio. 

Spiral structure of galaxies has eluded astronomers al- 
ready for a long time. Is it a quasi stationary density wave 
(Lin & Shu 1964, 1966)? Or is it a temporal phenomenon 
excited by internal or external disturbances and then swing 
amplified (Goldreich & Lynden-Bell 1965; Toomre 1981)? 
Or is it some kind of combination, like swing amplified den- 
sity waves? Until now theory and observations are inconclu- 
sive. To compound matters, there is the "anti spiral theo- 
rem" (Lynden-Bell & Ostriker 1967), which states that spiral 
structure cannot develop in a system governed by a time re- 
versible description. Consequently spiral arms cannot exist 
in an isolated pure stellar galaxy. Some kind of non time 
reversible action is needed, like external forcing or gas dis- 
sipation. The need for the latter mechanism is confirmed 
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observationally; in general when there is no gas there are 
no spiral arms. When studying spiral structure one thus has 
to include a gaseous component in the calculations. Doing 
that, calculations in three dimensions are required because 
for a typical disc the scaleheight of the stellar component is 
five times as large as the scaleheight of the gas layer. 

Some galaxies have bars, others don't. Ostriker & Pee- 
bles (1973) showed that a stellar disc can be stabilized 
against bar formation by adding a spherical (dark halo) po- 
tential. This has lead to the common believe that barred 
galaxies have a larger disc to halo mass ratio compared to 
non-barred galaxies, although there is no observational evi- 
dence to support this believe. Other theoretical and numeri- 
cal calculations (Efstathiou et al. 1982; Bottema & Gerritsen 
1997) investigating this matter reconfirm that a relatively 
larger disc mass results in a galaxy which is more unstable 
to bars. Therefore this ratio is inevitably an important pa- 
rameter in determining the appearance of a galaxy. Another 
effective parameter has to be the stellar velocity dispersion 
of a disc. Larger dispersions, expressed in a larger Toomre's 
(1964) Q value, lead to a more stable disc (Sellwood & Carl- 
berg 1984). In addition a substantial bulge makes a disc more 
stable against bars (Sellwood & Evans 2001). 

A proper investigation of disc galaxies involves a.o. dis- 
sipational and non linear processes. Therefore this subject 
is not easily handled by analytical calculations and one has 
to resort to numerical simulations. Results of such simu- 
lations can then guide the way to more detailed analyses 
of some specific aspects. Concerning numerical simulations, 
the combination of the tree code (Barnes & Hut 1986) for 
the collisionless part of a galaxy with the Smoothed Par- 
ticle Hydrodynamics (SPH) algorithm (Lucy 1977) for the 
collisional part allows the study of a galaxy containing gas. 
However, this combined TREESPH code (Hernquist & Katz 
1989) has rarely been used to study the detailed evolution 
and stability of isolated galaxies. It has mainly been applied 
to investigate interacting galaxies (Barnes & Hernquist 1991; 
Mihos & Hernquist 1994; Springel 2000), galaxy formation 
(Katz & Gunn 1991; Katz 1992) or barred galaxies specifi- 
cally (Friedli & Benz 1993, 1995). A preliminary numerical 
study of spiral structure has been performed by Elmegreen 
& Thomasson (1993) while Gerritsen (1997) and Gerritsen 
& Icke (1997, hereafter: GI97) made a detailed numerical 
analysis of the distribution and interaction of the stellar and 
gaseous components in an isolated galaxy. 

For a galaxy simulation one needs at least stars and gas. 
An optional bulge can be added and in order to make the ob- 
served flat rotation curves, a dark halo should be included. 
When gas is present, that cools and dissipates so that a star 
formation (SF) prescription is needed. That leads to gas 
depletion and to the formation of young stars with low ve- 
locity dispersions. SF also generates heating and stirring of 
the gas, preventing it from further collapse and so establish- 
ing a feedback loop which creates some kind of interstellar 
gas medium in equilibrium. These processes are too com- 
plex to be simulated from basic principles so that one needs 
a global description for the entities and processes in galaxy 
simulations. 

Concerning SF there is a clear observational result, 
namely the Schmidt (1959) law which states that the SF 
rate (SFR) is proportional to the gas density to a power n: 
SFR oc p n . Detailed observations and analyses by Kenni- 



cutt (1983, 1989, 1998) and by Ryder & Dopita (1994) show 
that this Schmidt law holds with n ~ 1.3 for a wide range 
of (surface) densities in galactic discs. Consequently this ob- 
servation and proper exponent has to be the outcome of any 
numerical SF prescription. The interesting result of GI97 is 
that they reproduce the right Schmidt law by using a star 
formation law based on a Jeans criterion. They do not im- 
pose the Schmidt law a priori, but it is the result of a more 
fundamental process. Their scheme is therefore expected to 
hold for a broad range of physical conditions. 

It is well known that in a real galaxy young stars and 
SF regions exert a definite action on the ambient ISM. The 
fierce radiation heats the gas and supernova (SN) action 
stirs up the ensemble of gas clouds. This SF feedback main- 
tains an equilibrium state of the ISM. When simulations lack 
such a feedback a catastrophic gas cooling occurs, forming 
super giant cloud complexes (Shlosman & Noguchi 1993). 
Therefore in the numerical calculations also a SF feedback 
prescription is needed. 

As mentioned above, Elmegreen & Thomasson (1993) 
made a numerical study of isolated galaxies to specifically 
investigate spiral structure, more or less as a function of disc 
to halo mass ratio. They find for a "rather massive" disc a 
grand design m = 2 structure which lasts for several rota- 
tion periods. For substantially less massive discs the spiral 
structure becomes flocculent. Though their investigation is 
very interesting it contains several caveats. The main one is 
that the development of a bar is prevented in at least two 
ways. Firstly by using an artificial Q barrier in the inner 
regions. Secondly the disc is stabilized by a large softening 
length (Romeo 1994) which has to be used to perform the 
two dimensional simulations. Therefore it is expected that 
the results of Elmegreen & Thomasson only partially ap- 
ply to a real galaxy. Results in the present paper differ from 
the study of Elmegreen & Thomasson mainly relating to the 
development of a bar in massive discs. 

My original intention was to investigate what happens 
to the morphology of a normal galaxy when the ratio of 
disc to dark matter changes. Inevitably, during the research 
other and additional problems emerged which necessitated 
a somewhat broader scope. Nevertheless this paper still fo- 
cuses on the global properties and global morphology of a 
galaxy and does not aim to go into specific details. Besides 
the abovementioned mass ratio another defining parameter 
for a galactic disc is the value of the stellar velocity disper- 
sion, or more or less equivalently, the thickness of the disc. 
In addition effects of the amount of SN feedback and amount 
of available gas have been considered. The investigations are 
restricted to late type galaxies for which NGC 628 serves as 
template. That galaxy is nearly face-on, has only a minor 
bulge, and has a well defined, non barred spiral structure. 

Section 2 of this paper describes the code and all the 
relevant prescriptions. The employed numerical scheme is 
based on the works of GI97, Gerritsen (1997), and Hern- 
quist & Katz (1989). Besides the nice asset of reproducing 
the Schmidt law, the scheme creates a multiple phase in- 
terstellar medium, equal to the situation in a real galaxy. 
GI97 performed their simulations for galaxies with low SF 
rates and then an equilibrium gas state can be achieved by 
UV heating only. When the SFR crosses a certain threshold 
it is necessary to include mechanical feedback in addition. 
This supernova feedback will be desribed in sufficient de- 
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tail. In Sect. 3 the model galaxy is described and in Sect. 4 
several properties and adaptations are discussed. To create 
a realistic disc in star formation equilibrium in Sect. 5 the 
necessary amount of SN feedback and external gas supply 
are determined. Section 6 describes what happens when the 
main parameters; disc to halo mass ratio and disc thickness 
are varied. When does a bar form and what is the global 
morphology? In Sect. 7 an analysis and discussion is given 
of the gas properties, like temperature distribution and gas 
velocity dispersion. Section 8 deals in some detail with the 
spiral structure and compares results from the simulations 
with observations. In Sect. 9 the theory of swing amplifica- 
tion is applied to the galaxies under investigation and finally 
in Sect. 10 a general discussion is given and a point by point 
summary of the conclusions. 



2 THE CODE 

2.1 Numerical technique 

The evolution of a galaxy is simulated using a hybrid N- 
body /hydrodynamics code (TREESPH; Hernquist & Katz 
1989). A tree algorithm (Barnes & Hut 1986; Hernquist 
1987) determines the gravitational forces on both the stel- 
lar and gaseous components of the galaxy. The hydrody- 
namic properties of the gas are modelled using the method 
of smoothed particle hydrodynamics (Lucy 1977; Gingold 
& Monaghan 1977). The gas evolves according to hydrody- 
namic conservation laws, including an artificial viscosity for 
the treatment of shocks. Each gas particle is assigned an 
individual smoothing length, hi, which determines the local 
resolution. Gas properties are found by smoothing over N 
neighbours within 2h, where N is typically 64 for a number 
of gas particles larger than 12.000. The acceleration (dl) of 
particle i is given by 



oi = VPi + cn viec -V*i, 

Pi 



(1) 



with P being the pressure, d*; vlsc the artificial viscosity term 
(Monaghan 1989), and $ the gravitational potential. For the 
evolution of the specific thermal energy (u) an expression is 
chosen as 



dm 11 



r- a 



(2) 



where (F — A) / p accounts for the non-adiabatic heating and 
cooling terms. To close the system of equations describing 
the evolution of the fluid an equation of state 



P = (7 - l)pu, 



(3) 



is adopted with 7 = 5/3 for an ideal gas. 

Time steps of the stellar particles were equal to 1.5 10 6 
years. For gas (SPH) particles a variable time step was used 
determined by the Courant-Friedrichs-Levy condition, with 
a maximum of 1.5 10 6 years. The softening of stellar parti- 
cles was not implemented by the usual softened force law: 
a oc (r 2 + e 2 ) -1 but by a smoothing kernel surrounding the 
stellar particle (Hernquist & Katz 1989) . In that way a bet- 
ter approximation is achieved of the real force law. 



2.2 The cooling 

The cooling term A in Eq. (2) describes the radiative cool- 
ing of the gas. For the simulations a constant composition 
of H I plus Helium gas is adopted. The "standard" cooling 
function of Dalgarno & McCray (1972) was used for this gas, 
parameterized as 
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where T is the temperature and nu the hydrogen density 
in cm The second term on the right hand side of Eq. 
(4) determines the cooling below 10 4 K, which is strongly 
dependent on the amount of ionization of the Hi gas. For 
an ionization parameter n e /riH of 0.1 (Cox 1990) a good 
approximation of Dalgarno & McCray's cooling function is 
found for a = 3.24, b = 0.085, and c = 3.0. 

The star formation rate is dependent on the amount of 
cooling below 10 4 K and the parameter a mainly determines 
this amount. For a lower ionization or for lower abundances 
the cooling will be smaller and the SFR lower (Gerritsen & 
de Blok 1999). In the calculations the cooling rate for all gas 
is set a priori equal to approximately the solar neighbour- 
hood value. The effects of metallicity changes in a galaxy 
have not been considered. 



2.3 Heating 

The heating of the gas (r) is described by two terms: 

r 



[erg s 1 cm -3 ] 



10 eGqUh + CcrUH, 



(6) 



where the first term represents photoelectric heating of small 
grains and PAHs (de Jong 1977; Wolfire et al. 1995) and 
dominates in the stellar disc region. In Eq. (6) e is the heat- 
ing efficiency and Go is the incident far-ultraviolet field (912 
to 2100 A) normalized to Habing's (1968) estimate of the lo- 
cal interstellar value (= 1.6 10~ 3 erg cm" 2 s _1 ). This heat- 
ing efficiency e w 0.05 for temperatures below 10.000 K. The 
second term in Eq. (6) represents the so called "cosmic ray" 
heating (Black 1987; Spitzer 1978). A value of 7.8 10~ 27 
erg s -1 has been used for the cosmic ray heating constant 
C cr . In a normal stellar environment cosmic ray heating is 
negligible compared to the grain heating. Yet in the outer 
regions of a galaxy, beyond the stellar edge, it is still small, 
but it is the only heating source. This term has been in- 
cluded after the work of GI97 to prevent too much cooling 
and star formation in those outer regions. 

Since the gas is heated mainly by FUV photons, the 
stellar radiation field has to be calculated. In the simula- 
tions each stellar particle effectively represents a stellar as- 
sociation. Such an association is formed instantaneously and 
consequently each stellar particle has an age and a belonging 
FUV flux according to that age. This flux has been calcu- 
lated using the population synthesis models of Bruzual & 
Chariot (1993) for a certain specified IMF. Knowing the 
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FUV flux for every stellar particle, the radiation field is cal- 
culated by summing the individual fluxes corrected for geo- 
metric dilution. No extinction corrections have been applied. 
The cooling of the gas increases by more than a factor 10 3 
for temperatures rising above 10 4 K. That temperature is 
therefore effectively an upper limit, virtually independent of 
the heat input. 



2.4 The star formation criterion 

Star formation on scales < giant molecular clouds is a very 
complicated process. It involves matters like heating/cooling 
balance, collapse, fractionization and sub fractionization, 
the transition of Hi into H2, magnetic fields, winds, SN ex- 
plosions, shocks, ionization, desintegration, etc, etc. GMCs, 
where the star formation takes place, have typical masses 
of a few times 10 5 to 10 6 Mq, densities > 10 3 cm~ 3 and 
lifetimes of 10 7 to 10 8 years (Shu et al. 1987; Bodenheimer 
1992). Masses of gas particles in the simulation range typ- 
ically from 5 10 5 to 10 6 solar masses, being comparable to 
GMCs. A star formation recipe is thus a prescription of con- 
verting GMCs into stellar associations and does not consider 
all the processes mentioned above. 

The formation of stellar associations (= stellar parti- 
cles) proceeds nearly equal as described by GI97. For con- 
venience a summary of this star formation process will be 
given. In SPH simulations the gas particles sample proper- 
ties such as density and temperature. For each particle one 
can therefore calculate the Jeans mass Mj 



3/2 



(7) 



where p is the density, c 3 the sound speed and G the grav- 
itational constant. If this Jeans mass is smaller than a pre- 
scribed critical mass M cr i t , which is of the order of the mass 
of a GMC, then the gas particle is in a cold and dense en- 
vironment and is considered unstable to star formation. So, 
the condition Mj < M cr i t is a first criterion a gas particle 
has to obey to form a stellar particle. 

Once a region is unstable, it takes a collapse time to 
form a stellar cluster. As a second criterion it has been cho- 
sen that a gas particle has to remain unstable for a time 
span longer than a fraction f c of the free-fall time t//: 



£span ^* fc X tff 



fc 



(8) 



As soon as an SPH particle has fulfilled both conditions (7) 
and (8) a part of its mass (M s ) is converted into a stellar 
particle with mass M* 



M* = e sf x M g . 



(9) 



This continues untill the mass of a gas particle would drop 
below 0.1 times its original mass, then the whole remaining 
gas particle is turned into a new star. The star particle is 
given the same IMF as all other star particles, and age zero. 
New star particles are then included in the calculation of the 
radiation field and the ambient gas will be heated. In this 
way a self regulating star formation process is established. 

GI97 made a thorough investigation of the effects when 
a specified SF parameter is changed. In short, the cooling 
function essentially determines the SFR and a change of the 




10 7 2xl0 7 3xl0 7 

Time (years) 



4xl0 7 



Figure 1. Graphical representation of supernova (SN) feedback. 
The dashed line and left axis give the number of SNs as a function 
of time following Bruzual and Chariot (1993). In the code the me- 
chanical SN feedback is implemented by giving a certain fraction 
(SN efficiency) of the belonging SN energy as a constant thermal 
energy to one SN particle. This is illustrated by full drawn line 
and right axis for a SN action from t = t^ cg to t = t en< j for egN 
= 50%. 



collapse factor f c changes the ratio of warm to cold gas. 
Changing other parameters has a negligible or only a small 
effect on observable structures and quantities. In the present 
study the star formation parameters are essentially fixed at 
a certain value. The critical mass M cr i t is equal to the value 
of a GMC, the IMF is Salpeter (0.1 -> 125 Mq), the cooling 
function is solar neighbourhood and e s f = 0.5. The collapse 
factor was determined at 0.5 such that the gas is not too 
clumpy and the warm to cold gas ratio is reasonable (see 
Sect. 5.1 and 7.2). 



2.5 Supernova feedback 

Besides the FUV feedback on the gaseous medium there 
is also the mechanical feedback by supernovae and stellar 
winds. A method to implement that feedback has been de- 
veloped by Gerritsen (1997) which I have largely taken over. 
Since a description has not appeared in the refereed litera- 
ture and because a number of changes have been made to 
the method, a detailed description now follows. 

Winds and supernovae from OB associations are ca- 
pable of blowing cavities in the surrounding ISM. In the 
standard model for such a process the mechanical energy 
is injected into a hot bubble which sweeps up the external 
medium by a dense shell (McCray & Kafatos 1987). Not all 
the energy is converted into kinetic energy of the expanding 
shell. Silich et al. (1996) demonstrate that most of the en- 
ergy is simply radiated away and only 5 to 10% ends up in 
the expanding shell. An amount of ~ 10 51 ergs of energy is 
provided by each supernova so that every star more massive 
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than 8 Mq supplies that amount to the mechanical feed- 
back process. Additional energy is provided by winds from 
massive and medium size stars (Lamers & Leitherer 1993). 

The more classical implementation of mechanical feed- 
back is by "kicking" the surrounding ISM. The gas around 
a young star cluster is then given a velocity perturbation 
directed radially away from the cluster. A free parameter 
in this implementation is the fraction of the released en- 
ergy (/„) that is put into the velocity perturbation. The 
feedback may have profound effects on the gas depending 
largely on the value of f v . For example, Navarro & White 
(1993) favour values of /„ < 0.1 in their cosmological simu- 
lations. But Friedly & Benz (1995) take of the order of 0.01 
for barred galaxies and Mihos & Hernquist (1994, 1996) use 
f v ~ 10~ 4 for simulations of interacting galaxies. Clearly 
there is a lot of arbitraryness as to the choice of f v . 

At present a different approach to implement mechan- 
ical feedback has been chosen. Following an idea by Heller 
& Shlosman (1994) a fraction of the mechanical energy is 
given as thermal energy to one gas particle during a cer- 
tain amount of time. In this way a hot bubble is simulated 
by making one SPH particle hot; from now on called SN 
particle. It is further left to the TREESPH program to dis- 
tribute the energy to the surrounding gas particles like a 
pseudo-point explosion. 

For the implementation in the code it has to be realized 
that one is dealing with particles representing GMCs and 
whole star forming regions. When a gas particle (GMC) is 
deemed to form a star particle (SF region) half of its mass 
is transformed into a star particle with zero age and belong- 
ing IMF. The other half of the gas particle is designated as 
SN particle possessing the mechanical energy appropriate 
for the SF region. For a period of ~ 3 10 7 years SN activity 
is present in the SF region and for such a time the SN par- 
ticle is both, fixed to the parent star particle, and kept hot. 
To illustrate matters, in Fig. 1 the SN energy is given as a 
function of time according to the Bruzual & Chariot pop- 
ulation synthesis code. On top of that figure the simplified 
SN particle energy is displayed. For a small time span after 
star formation SN activity is not yet switched on and the 
SN particle is kept at a temperature of 10 4 K. After that a 
fraction esn of the total SN energy of the SF region is given 
as thermal energy to the SN particle. This total SN energy 
is calculated from the total mass of the new star cluster, the 
assumed IMF, and the fact that all stars go SN for a mass 
> 8 Mq. After a time tfjf no more energy is supplied to 
the SN particle and it is allowed to cool again which it will 
do quickly in practice. 

In reality there is a time variability of the SN action 
(Fig. 1) and a poorly known energy contribution from winds. 
These details and uncertainties have been included in the 
free parameter £sn- After the SN phase the gas particle re- 
turns to being a normal SPH particle which may undergo a 
next event of star formation if the conditions are right. Most 
simulations in this paper are performed for esn = 25%, and 
for a regular SF where half the gas mass is transformed into 
stars, the belonging temperature of the SN particle amounts 
to 1.25 10 7 K. In this, albeit artificial, way a three phase 
ISM is created containing cold, warm, and hot gas (McKee 
& Ostriker 1977). 

Gerritsen (1997) did some testing of this method on iso- 
lated galaxies. It appeared that realistic cavities are created 
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Figure 2. The total gas surface density of NGC 628 which has 
been used in the model galaxy. H I and H2 have been added and 
multiplied by 1.4 to account for Helium. 



in the ISM, where the holes are smaller in the central denser 
regions and larger near the rim of the gas distribution. The 
size and the size distribution of the holes does not depend 
on the number of gas particles. This demonstrates the relia- 
bility of the method and shows it can be used for a large gas 
density range even at a limited resolution. Compared to a 
situation with only FUV heating of the gas, adding the me- 
chanical feedback changes the gas appearance to a different 
and larger design. 



3 THE GALAXY 

All the calculations have been performed for one model 
galaxy with different internal compositions. This galaxy con- 
sists of four components: a stellar disc, a gas layer, a bulge, 
and a dark halo. This disc and gas layer are made up of 
stellar and gas particles while the bulge and halo are imple- 
mented as rigid, "dead", potentials. 

An exponential stellar disc is assumed with a locally 
isothermal z-distribution (van der Kruit & Searle 1981) 



Pd{R,z) 



pQe R ^ h sech~ 



(10) 



where h is the scalelength of the disc and zo the thickness 
parameter which was constant as a function of radius. If 
the disc were isolated the vertical stellar velocity dispersion 
(Vz) 1 ^ 2 depends on the stellar surface density <r» as 



2,1/2 



y/TrGz a*(R). 



(11) 



In reality adjustments are necessary to this relation because 
the other galactic components influence the stellar velocity 
dispersion. For the setup of the galaxy only the gas has been 
taken into account such that in Eq. (11) a t (R) is replaced 
by a*(R) + cTgas(-R). During the initial stages of the simula- 
tions the stellar disc still has to settle in the bulge and halo 
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Table 1. Galaxy parameters 



Disc 

V max 

(km s- 1 ) 


Disc 
mass 

(10 9 M ) 


Disc central 
surf. dens. 

(M 0P c- 2 ) 


{Vz)r^o for 
2 = 900 pc 
(km s- 1 ) 


Disc 
{Ml L B ) 


Halo 

Rcorc 

(kpc) 


Halo 

V max 

(km s- 1 ) 


Bulge 
mass 
(10 9 M ) 


Bulge central 
surf. dens. 

(M 0P c- 2 ) 


126 


45.5 


367 


67 


2.05 
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185 


2.4 


2562 


150 


64.5 


520 


80 


2.90 


5 


180 


3.4 


3630 


165 


78.0 


629 


88 


3.51 


11.2 


207 


4.2 


4390 


180 


92.9 


749 


96 


4.18 


23.4 


264 


5.0 


5230 



Distance = 11.5 Mpc 

Total light in B = 22.2 10 9 L Q (uncor.) 

Total gas mass = 18.8 10 9 M 

Disc scalelength = 4.5 kpc 

Bulge scalelength = 397 pc 

Disc mass / bulge mass = 18.7 



potential, but in general such a disc adaptation is small, cer- 
tainly if the disc is relatively massive. The dispersions in the 
radial and tangential directions are related to (t; 2 ) 1 ^ 2 as 



(v^=0.6g(R){v R )^, 



and 

/ 2\l/2 
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(13) 



Here g(R) is a function close to 1.0 and will be discussed 
in the next section. The tangential dispersion (iig,) 1 ^ 2 is 
related to the radial dispersion through O and k, which 
are the orbital and epicyclic frequencies. The average ro- 
tation of the stars (v„) is given as the testparticle rotation, 
vt = \J (d<&/dR)/R diminished with the asymmetric drift: 



2 2 / 2 \ 

v t - v, = (v R ) 



U* OR ) 



(14) 



using a plane parallel epicyclic approximation. At the radii 
where the stellar velocity dispersion becomes larger than the 
rotation, Eq. (14) does not hold any more. For these inner 
regions an interpolation by hand has been made to R = 0. 
Stellar particles are distributed according to Eq. (10) and 
are given a rotation according to Eq. (14). Dispersions in 
the R, z, and © directions follow from Eqs (11) to (13) and 
are drawn from Gaussian distributions. 

The gas layer has a radial functionality as observed for 
NGC 628. The thickness of the gas layer is determined by the 
surrounding potentials of the other components and the gas 
velocity dispersion. For the latter a constant, isotropic value 
of 8 km s _1 has been adopted. To calculate the thickness of 
the gas layer the (ad hoc) description of Bottema (1996) has 
been used. 

The bulge is spherical with a projected exponential ra- 
dial distribution, Ob = Co c~ r ^ hb resulting in a spatial dis- 
tribution of the bulge density p^ of 



Pb 



A. 



*(£)■ 



(15) 



Here hh is the scalelength of the bulge and Ko a modified 
Bessel function. A cutoff to the bulge density is taken at 
eight bulge scalelengths. For the dark halo the usual pseudo 
isothermal sphere is assumed with a density ph of 



Ph = Po 



1 + 



R 2 



and belonging rotation velocity v h of 



Rc 
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-arctan 
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where i? C orc is the core radius of the halo related to the 
maximum rotation wjj,. 



by 



^/I^GpJrTo 



(18) 



As template for the model galaxy the beautiful face-on 
spiral NGC 628 is used. It has a regular spiral structure and 
certainly not a bar. This gives the possibility to compare 
the structure resulting from the simulations with an actual 
galaxy in an unambiguous way. It is not intended to make 
a model of NGC 628 in all its details, but only to generate 
a galaxy which has roughly the same dimensions. Photom- 
etry of NGC 628 is presented by Shostak & van der Kruit 
(1984) and by Natali et al. (1992). The first authors give a 
disc scalelength for the Kodak J emulsion of 85" while the 
second authors give hd = 73" in the R-band. As a compro- 
mise I shall use a value of 80" which gives for an adopted 
distance of 11.5 Mpc a value of h d = 4.5 kpc. A photometric 
disc profile with that scalelength has been subtracted from 
the photometry of Natali et al. and to the remaining central 
bulge light an exponential has been fitted. From that follows 
a bulge scalelength of 7" 12 = 397 pc and the surface bright- 
ness of bulge and disc are equal at a radius of 15" 3. This 
photometric analysis gives the absolute scales and bright- 
ness ratio of bulge and disc but not yet the masses of these 
components. 

NGC 628 is observed in Hi by Wevers et al. (1986) who 
give the radial H I gas density. This profile is supplemented 
with the H2 radial profile as observed by Combes & Bec- 
quaert (1997) and the total is multiplied with a factor 1.4 to 
account for Helium. The resulting radial gas profile as it has 
actually been used is given in Fig. 2. Note that although the 
addition of H2 in the centre looks quite impressive, it only 
contributes 8.5% to the total gas mass of 18.8 10 9 Mq. 

Because NGC 628 is face-on its rotation curve cannot 
be measured. Yet for its absolute luminosity and Hubble 
type the Tully-Fisher relation of Rubin et al. (1985) gives 
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Figure 3. For the photometry and gas content of NGC 628 a 
decomposition has been made of an adopted rotation curve reach- 
ing a flat level of 200 km s _1 (full drawn line). Four different 
disc masses are considered corresponding to maximum disc rota- 
tions of 126, 150, 165, and 180 km s _1 (long dashed lines). The 
combined halo plus bulge rotation (short dashed lines) is then ad- 
justed to fit the adopted total rotation. For each disc contribution 
the stellar rotation is given by the asterisks. 

a maximum rotation of approximately 200 km s _1 . This 
value has been adopted as the exact maximum rotation for 
the model galaxy together with the observational establish- 
ment that rotation curves are flat in the outer regions of a 
galaxy. For the exponential disc rotation curves have been 
calculated (Casertano 1983) scaled to four different maxi- 
mum rotational values of 126, 150, 165, and 180 km s _1 . 
These four stellar discs all have a certain mass and mass- 
to-light ratio. Combined with each disc is a bulge density 
distribution where equal M/L ratios of bulge and disc have 
been assumed. For the four cases a dark halo has been cre- 
ated such that the total rotation is approximately flat at 
200 km s _1 . The result of this fitting procedure is shown in 
Fig. 3 and numerical values are presented in Table 1. Also 
given in Fig. 3 is the actual stellar rotation which results 
after the subtraction of the asymmetric drift, in this case 
for a zo value of 900 pc. 



4 PRACTICAL PROPERTIES AND 
ADAPTATIONS 

4.1 Main parameters 

To gain insight into the possible morphologies of a disc, 
the main structural parameter to be investigated is the ra- 
tio of maximum disc rotation to total rotation: «S??/ w Sbs X - 
Though it has already been known since the calculations 
by Ostriker & Peebles (1973) that a galaxy where the disc 
approaches its maximum disc limit is unstable to bar for- 
mation, (see also Bottema & Gerritsen 1997) that notion 
should be confirmed for an actual 3D disc with gas and 



star formation. Inevitably related to the matter of stability 
must be the thickness of the disc. A thicker disc has a larger 
vertical velocity dispersion. The dispersions in the three di- 
rections are related, which means that a thicker disc has a 
larger Q value and should at least be more stable to local 
distortions. Therefore zo is the second main parameter to be 
investigated. 

The precise amount and the process by which mechan- 
ical SN energy is transferred to the gas is not known. Con- 
sequently the influence of the parameter esn will be con- 
sidered in some detail. A fourth parameter has been added 
during the course of the simulations. For a resulting SFR 
of 4 MQyr -1 it appeared that the gas reservoir becomes 
depleted pretty fast. To create a constant situation over a 
reasonable time span, gas has been added at several rates. 
Summarizing, the main parameters investigated in this pa- 
per are: z o, £sn, and gas supply. 

4.2 Particle masses and numbers 

For most simulations the number of stellar particles was 
150.000 and the number of gas particles 16.000, though some 
simulations are presented using 40.000 gas particles. In Ta- 
ble 1 the total stellar disc masses are given for the various 
maximum disc rotations. From = 126 km s _1 to 180 

km s _1 the stellar particle masses then range from 0.3 10 6 
to 0.6 10 6 M . For a total gas mass of 18.8 10 9 M Q and 
16.000 gas particles, the mass of one gas particle amounts to 
1.18 10 6 Mq. During the star formation process new stel- 
lar particles are created with masses of | , \ , and | of the 
mass of a gas particle, being comparable to the masses of 
the already existent stellar particles. 

4.3 The star formation rate 

The SFR as defined here is the total amount of gas mass 
which is transformed into stellar mass, per time unit. In the 
code no stellar mass is returned to the gas reservoir. For an 
actual galaxy with a standard IMF, approximately 70% of 
the gas that has formed stars is eventually locked up in old 
stars and stellar remnants. The other 30% is returned, of 
which most in a short time scale after SF. Because this time 
scale is so short the non returning of gas in the numerical 
simulations is a good approximation. However, the star for- 
mation rate that could be observed is then not equal to the 
SFR as defined above. For a lock up rate of 70% the "ob- 
servable SFR" is defined as 1/0.7 times the SFR and can as 
such be compared with observations. Kennicutt (1983) has 
determined an "observable SFR" of 3.8 MQyr" 1 for NGC 
628, when it is put at a distance of 11.5 Mpc. There is a 
systematic error of approximately a factor of two associated 
with this value of 3.8 MQyr -1 mainly caused by the uncer- 
tainty of the IMF. 

4.4 Stochastic heating 

An average galactic disc is essentially collisionless. In a sim- 
ulation a disc is made up of far less particles than the num- 
ber of stars in a real galaxy. As a result encounters between 
particles will occur in simulations and the particle ensem- 
ble will be heated up. This process of particle scattering is 
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Figure 4. The ratio of vertical to radial stellar velocity disper- 
sion of the model galaxy. The dynamical equilibrium is indicated 
for two cases: -u™™ = 150 km s _1 with z = 900 pc and d™^ = 
126 km s _1 with zn = 675 pc. Simulations are started with the 
functionality given by the "input" line with a z/R dispersion at 
R = of 0.9, 0.8, 0.8, and 0.7 for v%?£ = 180, 165, 150, and 126 
km s _1 respectively. 



commonly referred to as stochastic heating. In practice this 
means that a model disc made up of a limited number of 
particles will heat up, get thicker, and become more stable 
during a simulation. 

In a real galaxy which contains gas, stellar heating may 
be generated by two processes. At first heating by encoun- 
ters of stars with molecular clouds (Spitzer & Schwarzschild 
1953; Villumsen 1985). Secondly a more effective mechanism 
seems to be heating by transient spiral structures (Barbanis 
& Woltjer 1967; Carlberg & Sellwood 1985) though in prac- 
tice a combination of both mechanisms is likely taking place. 
In a simulation with a limited number of stars, and with spi- 
ral arms and cloud complexes both, stochastic heating and 
real heating will take place. The problem is to disentangle 
the two. 

This disentanglement can be done by considering the 
change of the stellar velocity dispersion. From basic princi- 
ples one can derive for stochastic heating that 

such that the change of the square of the velocity disper- 
sion is proportional to a function of the number of particles 
N. During a simulation the square of the velocity dispersion 
will increase by stochastic heating and by heating of clouds 
and arms. If the same simulation is repeated for different 
numbers of stellar particles both heating processes can be 
disentangled. Another way to study the matter is by com- 
paring the same simulations with and without gas, as long 
as the gasless simulation produces a smooth featureless disc. 
Using these methods, a small investigation was done for a 
simulation using = 150 km s _1 and zo = 675 pc. It 

appears that for a total of 180.000 stellar particles the disc 
does not show stochastic heating at a radius of one scale- 



length and shows it barely at two scalelengths. As a final 
compromise a number of 150.000 stellar particles has been 
used for all discs. This means that beyond a radius of ap- 
proximately two scalelengths a small amount of heating will 
be present, not occurring in real galaxies. 



4.5 The softening length 

When performing 2D simulations of discs, the thickness is 
generally artificially implemented by taking a stellar soften- 
ing of size comparable to the scaleheight of the disc. The 
particles then do not meet in the plane but a disadvantage 
of this adaptation is a considerable additional stabilization 
compared to a 3D situation with the same dimension and 
kinematics (Romeo 1994). 

For 3D simulations considerations on the size of the soft- 
ening length have been made by Hernquist (1987). It appears 
that in order to achieve a good approximation of the force 
computation the softening length should be smaller then 
A, being the mean interparticle distance at the half mass 
radius. For an exponential disc Hernquist gives an approxi- 
mation of A w 2.63(iV/fo 2 zn)~ 1/ ' 3 which for the present case 
of NGC 628 with N = 150.000 gives A « 120 pc. One may 
consider other evaluations of A, for instance the mean inter- 
particle distance at the half mass height above the plane. In 
that case A ranges from ~ 110 pc at the centre to ~ 230 pc 
at R — 2h. When a bar develops, however, A may become 
even smaller in the bar region. Therefore, for the galaxy un- 
der investigation in this study the softening length should 
be smaller than approximately 100 pc. 

To get a lower limit on the softening length is more 
problematic. If e is too small, strong interactions between 
particles will occur. That is especially problematic when 
particles with different mass are present; the lightest particle 
will then be ejected, carrying with it large amounts of energy 
and angular momentum. In addition, for close particle inter- 
actions the time step may be too small to achieve a correct 
orbit integration. To investigate matters, pure stellar simu- 
lations have been performed for a case where a moderate bar 
develops (Fig. 10, bottom left) with e ranging from 2 to 200 
pc. For e between 10 and 200 pc the resulting structure and 
stellar kinematics is practically equal. For e < 10 pc particles 
get ejected from the galaxy (particles have masses differing 
~ a factor two) which changes the eventual morphology. It 
is therefore concluded that for the present calculations there 
is a safe range between 10 to 100 pc for the stellar softening 
length. It has been fixed at 20 pc for all simulations. In no 
case there will then be artificial stabilization by a too large 
softening length. 



4.6 The quotient of vertical to radial stellar 
velocity dispersions 

In first instance this quotient was put at 0.6 for all radii as 
an extrapolation of the number for old disc stars in the solar 
neighbourhood. Consequently the function g(R) in Eq. (12) 
= 1.0 for all radii. During the initial stages of the simula- 
tion, for the more massive discs, there appeared to be some 
outflow of stars from the central regions. A solution to this 
problem was found after comparing the final quotient for 
some realistic simulations with the original input (see Fig. 
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Figure 5. The resulting galaxy after 1.5 Gyrs of simulation. In this and following similar figures the left column gives the simulated 
V-band image, the middle column gives the stellar mass image and the right column the gas mass. The disc has a maximum rotation 
of 126 km s _1 and zq of 675 pc while gas was supplied at a rate of 4 MQyr -1 . From top to bottom row the mechanical SN feedback 
efficiency was at 0%, 5%, 25%, and 50%. In case of none or low SN feedback unstable giant SF regions appear. More SN feedback results 
in a galaxy which is in SF equilibrium having an appearance not unlike that of a real galaxy. High SN feedback destroys the coherent 
gas spiral structure. 
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4). At the end of the simulations the velocity dispersion el- 
lipsoid is always more spherical in the inner regions, while 
it remains at ~ 0.6 near two scalelengths. The initial model 
galaxy has been adapted to have a function g(R) slightly 
variable as a function of radius. As shown in Fig. 4, the cen- 
tral quotient was put at 0.7 (for i>™ = 126 km s _1 ), at 0.8 
(for = 150 and 165 km s" 1 ) or 0.9 (for = 180 

km s _1 ) and decreased linearly to 0.6 at a radius of 10 kpc. 
For such a setup the initial central instability has completely 
disappeared. For a disc with maximum disc rotation of 126 
km s _1 or lower this adaptation was, in principle, not neces- 
sary; the galaxy slowly changes to the final quotient without 
noticeable redistribution of matter. 



5 CREATING A GALACTIC DISC IN 
EQUILIBRIUM 

5.1 The amount of SN action 

The first parameter of which the effect has been investigated 
is the SN efficiency. To that aim a realization of the model 
galaxy has been chosen which is close to what observations 
of galactic discs suggest. Therefore a thickness zq of 675 pc is 
assumed (van der Kruit & Searle 1982) and a disc maximum 
rotational contribution of 63% (Bottema 1993). By no means 
the statement is given that this galactic constitution is the 
closest to reality; only that it is a reasonable assumption 
and suitable to test the effect of the amount of SN feedback. 
During the simulations gas has been supplied at a rate of 

4 MQyr -1 being approximately equal to the amount of gas 
depletion by SF. 

Results of esn = 0%, 5%, 25%, and 50% are given in Fig. 

5 where for each case the simulated V-band image, stellar 
mass image, and gas mass image is presented. In fact, Fig. 
5 is pretty self explainable. If the SN feedback is absent 
or small, dense clumps of gas and young stars are created, 
which dominate the V-band and gas image. On these clumps 
cold gas continues to accumulate and because the density is 
so large, the cooling outruns the heating by the FUV field of 
young stars. It should be kept in mind that cooling proceeds 
proportional to the square of the gas density (Eq. 4) and 
heating proportional to the first power of the density (Eq. 6) 
and a runaway process as observed might indeed be expected 
to occur. For an esn of 25% the clumping process has almost 
completely disappeared and a number of spiral arms appear 
in the stellar light, mass, and gas distribution. Certainly in 
the gas structure these arms have a certain width contrary 
to the lower esn simulations where the gas arms are thin 
and have a threaded appearance. Increasing the supernova 
feedback to the high value of 50% results in the creation 
of multiple holes in the gas structure. The feedback is now 
so violent that spiral arms in the gas image have all but 
disappeared. As a result also the strength of the spiral arms 
in the light and stellar mass images decreases. 

For the different SN feedback situations the resulting 
SFR as a function of time is presented in Fig. 6. In case of 
low £sn and the generated severe clumping the SFR keeps on 
increasing, evidencing the runaway star formation process. 
If this would continue, the galaxy gets quickly depleted of 
its gas and the SF activity would subside again. A larger 
SN feedback results in a lower SFR, though above esn = 
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Figure 6. The star formation rate during the simulations pre- 
sented in Fig. 5. For low efficiencies of the SN feedback (0 and 
5%) SF is boosted in numerous giant SF regions. A larger SN 
feedback results in an equilibrium situation at which the SFR 
reaches a constant level determined by the gas supply and cool- 
ing function. 

25% the SFR remains more or less constant at a level of 4.2 
MQyr -1 . Corrected for stellar mass loss (Sect. 4.3) this gives 
6 MQyr -1 which is in good agreement with the observed 
value for NGC 628 (Kennicutt 1983). There is no evidence 
for the process of stochastic self propagating SF by super- 
novae explosions (Mueller & Arnett 1976; Elmegreen 1992). 
On the contrary, more SN action decreases the amount of 
SF when the SN action is low, while for higher regimes of 
SN action the SFR seems to be almost independent of the 
amount of SN feedback. During the simulations star forma- 
tion occurs mainly by cooling and dissipation of the gas. SN 
action prevents a runaway SF process and influences the gas 
morphology. 

An investigation has been made of how the structures 
and relations of Fig. 5 depend on other parameters. It does 
not depend on the value of M cr i t (Eq. 7). There is a depen- 
dence, however, on the value of the collapse factor f c (Eq. 8). 
If the time span before star formation is increased from the 
nominal value of 0.5 x tff the relative amount of cold (T < 
500 K) gas increases and the clumping process sets in at 
higher SN feedback values, and vice versa. Since the ratio of 
cold to warm H I gas is as yet undetermined it is impossible 
to pinpoint down the exact value of esn (see Sect. 7 for an 
ample discussion) . Still, following indications of the observa- 
tions leading to the adopted values for the parameters, the 
relation between the morphology in Fig. 5 and the value of 
esn is likely accurate to within a factor of two. 

5.2 Different amounts of gas supply 

As noted before, it became evident during this investigation 
that when having a gas consumption rate of ~ 4 Mgyr -1 
for a total gas amount of 18.8 10 9 Mq, the gas gets depleted 
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Figure 7. The star formation rate during simulations with dif- 
ferent amounts of gas supply. Eventually the SFR will adjust itself 
to the amount of gas that is being added. 



pretty fast. In the mean time, it takes a while, lets say some 
0.8 10 9 yrs, before a stable galaxy emerges (see Fig. 6). To 
have a constant gas reservoir over at least this adjustment 
time it was decided to create the possibility to add gas at 
a certain rate. For a real galaxy it is also expected that 
gas gets accreted, although the way how is (still) unclear. 
Therefore this gas suppletion is less artificial than it might 
seem. 

Care has been taken that the addition of gas does not 
alter the available gas structure. For example, gas should not 
be supplied at those positions where it has just recently been 
removed by star formation. Gas particles have been added 
randomly by selecting another gas particle and placing the 
new gas close to that selected particle. From this selection 
are excluded gas particles with temperatures above 20.000 K 
or below 5000 K. Also excluded are, of course, SN particles 
or particles already Jeans unstable. It became clear that 
when gas is added at places beyond the stellar disc, that 
only increases the amount of gas over there, taking more 
computing time, and not replenishing the gas where it is de- 
pleted. Therefore no gas has been added at radii larger than 
25 kpc. It was originally expected that extra gas at large 
radii would slowly migrate inwards by dissipation and vis- 
cous interaction. That process, however, did not take place. 

The morphology and SFR of a galactic disc has been 
investigated in cases of gas addition of 0, 4, and 8 MQyr -1 
during a period of 1.5 10 9 years. A disc is chosen equal to the 
one where the amount of SN action was investigated, with 
a fixed £sn of 25%. A SFR at a constant level is achieved 
when ~ 4 MQyr -1 is added. When no gas is added the 
SFR approaches 4 MQyr -1 after 3 to 4 10 8 years from the 
start of the simulation, and declines gradually thereafter (see 
Fig. 7). At the end of the simulation the SFR has decreased 
to 2 MQyr -1 . As can be seen in Fig. 7, when gas is added at 
8 Mq yr -1 the SFR keeps on increasing to reach eventually 
the level of gas input. In this way one can draw an impor- 



tant though not unexpected conclusion: the SFR in a disc 
adjusts itself to the amount of gas that is being supplied. 
The time scale on which this happens is for the present sim- 
ulation ~ 1 Gyr, but will likely depend on the numerical 
resolution. When the images of the galaxy are inspected at 
1.5 Gyrs after the start of the simulation it appears that the 
morphology is only slightly dependent on the amount of gas 
supply and consequently also only slightly dependent on the 
actual SFR. More gas and higher SFR results in a somewhat 
stronger and less regular spiral pattern. When 8 MQyr -1 is 
added also a few dense gas clouds develop. 

5.3 Comparison with observations 

The V-band images as presented in the left column of Fig. 5 
are made to resemble a galaxy as depicted in an optical 
photograph. To that aim the intensity is given on a square- 
root gray scale and cut off (saturated) above a certain level. 
There is always some subjectiveness involved in selecting the 
scale and cutoff and for the best result I aimed at achieving 
the best overall contrast. In reality, however, this will be 
equal to the process of developing and printing of actual 
photographs. A shortcoming of the simulations is the lack 
of dust. Certainly the young SF regions in real galaxies are 
often (at least partly) enshrouded in some amount of dust, 
a situation that cannot be reproduced. Another problem is 
that every new star particle (= SF region) has an IMF going 
up to stars of 125 Mq which will completely dominate the 
light in every passband. In a real galaxy there are only a few 
SF regions which contain such massive stars. To compensate, 
at least in part, for the lack of dust and for the discrete 
appearance of very massive stars, the V-band image is made 
for all the stars at an age 3 10 7 years older than the actual 
age of the stars. The net effect of this procedure is an image 
less dominated by young star clusters. Admittedly, there is 
some artificiality in creating a simulated V-band image to 
compare with photographs, but there is no way around this. 

When there is no SN action the V-band image is com- 
pletely dominated by bright Hn regions and to my knowl- 
edge such a galaxy does not exist. For esn = 5% the optical 
image is also very knotty and is barely consistent with the 
image of NGC 628 on, for example, page 99 of the Shapley- 
Ames catalog. When £sn is larger the optical images are not 
unlike images of real galaxies. The problem with the stellar 
mass images is that these cannot be observed. Even in the 
near infrared passbands young populations will supply ex- 
tra light, e.g. have a smaller M/L ratio compared to an old 
population. 

Contrary to the optical and stellar mass image, the im- 
age of the gas can actually be compared with observed H I 
structures. High resolution Hi observations of nearby late 
type galaxies are presented by a.o. Braun (1995) and Kam- 
phuis (1993). Such a high resolution is necessary because 
one would like to compare the small scale gas structure and 
in particular the size and distribution of holes created by SN 
feedback. For example, Fig. 5 lower right panel shows lots of 
small holes and cavities in the inner region, all created by SN 
explosions. In addition, a particular feature of SN feedback 
are holes blown into the ridge of spiral arms. This can be 
witnessed in the same image at for example positions (x,y) 
= (5,-15) and (-4,19). The observations of the gas structure 
do never show the knotty and threaded gas structures as for 



12 Roelof Bottema 



£sn = and 5%. The best resemblance between simulations 
and observations is for esn around 25%, though also more 
amorphous gas structures resembling the structure for esn 
= 50% are observed. A complicating factor is that some of 
the nicest H I spiral structures are observed in galaxies with 
grand design spiral arms obviously caused by an interaction 
with other galaxies, for example M81 (= NGC 3037). In 
such cases gas becomes aligned with the gravitational spiral 
disturbance, influencing its global appearance. This should 
be taken into account when comparing the simulated gas 
structures with real observations. 



5.4 Whence spiral arms? 

The spiral structure which is generated in the present simu- 
lations of an isolated galaxy is the result of an equilibrium of 
four actions. On the one side cooling and dissipation accu- 
mulates the gas in threads and clumps. This is counteracted 
by FUV heating and SN feedback. If the latter action is too 
small the gas collapses completely into a few giant SF re- 
gions. On the other hand, if the action of heating and SNs 
is too large, warm gas get dispersed over the disc prevent- 
ing the formation of any coherent structure. Except for this 
amount of mechanical feedback the resulting gas and mass 
structures do not depend on details (like IMF, M cr i t , etc) 
of the star formation process (see Sect. 5.1 and GI97). It is 
therefore likely that these structures do not depend on the 
exact star formation recipe. 

The behaviour of the gas is the main cause of structure 
formation. As a result the underlying stellar mass distribu- 
tion becomes disturbed and also in the stellar disc spiral 
arms are formed. In this way gas mass enhancements more 
or less swing amplify the stellar mass into a spiral structure. 
But this is probably not the whole story. If it were only the 
need of mass accumulations to generate spirals, then it is 
expected that the low esn situation with its specific clumps 
would generate the strongest spirals. Obviously that is not 
the case. The gas enhancements really have to be wound up 
in spirals by the differential rotation to create an accompa- 
nying mass enhancement. A situation is established where 
shearing gas filaments swing amplify the underlying mass 
distribution into a transient small scale density wave. The 
latter appears as a spiral arm. More insight is gained when 
the simulation is played as a movie on the TV screen. It can 
then be noticed that coherent spiral features develop and 
exist for approximately half a rotation period. Spiral arms 
are "wound" into existence and disappear by winding and 
stretching too much or by interaction with other spiral fea- 
tures. Such interactions also create the typical bifurcation 
patterns of the arms often seen in real galaxies. Not created 
in the present simulations are the grand design m=2 spirals 
witnessed in a few exemplary galaxies. These matters will 
be further investigated and discussed in sections 8 and 9. 



6 THE BAR INSTABILITY 

6.1 The appearance as a function of disc mass 

An amount of SN feedback has now been determined which 
results in a good resemblance with observed structures in 



a disc while for a gas addition of 4 MQyr~ the gas reser- 
voir remains constant. Fixing these parameters and fixing 
the thickness parameter zo at 900 pc the resulting morphol- 
ogy for different disc contributions to the total rotation has 
been investigated. Images of the galaxy after 1.5 10 9 years of 
simulation are presented in Fig. 8 for t^Jc °f 126, 150, 165, 
and 180 km s _1 . As for Fig. 5, Fig. 8 speaks pretty much for 
itself. Discs with a moderate rotational contribution of 126 
and 150 km s _1 evolve into a regular spiral galaxy. For v™^* 
= 150 km s _1 there appear to be less arms which are more 
pronounced than in the case of «3J?<* — 126 km s _1 . However, 
the images presented in Fig. 8 are only one snapshot at one 
particular time. During the simulations arms come and go 
and are variable in strength. Judging the whole simulation 
from t = to t = 1.5 10 9 years the vfitc = 150 km s" 1 
case has also on average a slightly more pronounced spiral 
structure. 

Discs with a larger rotational contribution obviously de- 
velop a bar. For v™^* = 165 km s _1 the disc starts out by 
forming a few large and broad spiral structures, a situation 
which continues to ~ 1.2 10 9 year. Then quite suddenly the 
bar appears. For v™^* = 180 km s _1 the simulation also 
sets out by forming broad dominant spiral features which 
already after ~ 0.5 10 9 years settle into a barred structure. 
When the bar is present gas is transported inwards rapidly 
along the bar. The same conclusion was reached by Athanas- 
soula (1992) for simulations of gas moving in a rigid barred 
potential. That study shows that the gas is transported in- 
wards along shocks appearing on the leading side of the bar. 
The present calculations give the same results, though the 
resolution in not good enough to investigate details of this 
process. 

In the bar the gas density is large, cooling is fast and 
new stellar particles form in appreciable quantities. So much 
young stars form at some positions and certainly at the cen- 
tre that FUV heating outweighs the cooling and gas remains 
at ~ 10 4 K while becoming ever more dense. When gas parti- 
cles accumulate within the softening length numerical prob- 
lems occur. To avoid these, the gas was forced to form new 
stars above a density limit of 10 Mq pc -3 even though the 
star formation criteria of subsection 2.4 were not yet satis- 
fied. It was not intended and it is not possible to simulate the 
processes going on at the very centre of the bar region. What 
the simulations do show is that gas is transported inwards 
along the bar, where it forms stars at a considerable rate. 
As a consequence a concentration of young stars emerges at 
the inner few kpc which takes on a more or less spherical 
shape. This can be called the development of a bulge. For 
the whole galaxy the radial luminosity profile changes from 
exponential to Freeman type II (Freeman 1970) where the 
shoulders of the profile coincide with the bar. Getting back 
to Fig. 8, one can notice the effect of the inwards gas trans- 
port. For = 180 km s _1 the bar has existed longer than 
for VjiJc = 165 km s _1 . As a consequence in case of v™^* 
= 180 km s^ 1 the region around the bar has already been 
severely depleted of gas. The results of the present simula- 
tions of the barred galaxies are identical to those of Friedli 
& Benz (1995). They also specifically report the gas inflow, 
SF along the bar and in the centre, and the formation of a 
bulge. 

But, do we need the gas to form or not to form bars? 
To investigate this, the simulations presented in Fig. 8 have 




Figure 8. For a fixed SN efficiency of 25%, gas suppletion rate of 4 Mqjt: -1 , and zq = 900 pc, the resulting galaxy after 1.5 Gyrs of 
simulation is presented for four values of the maximum rotation of the disc (126, 150, 165, and 180 km s -1 ). Less massive discs evolve 
into a regular spiral galaxy, while more massive discs develop a bar. In case of a bar, gas is transported inwards as can be noticed in the 
right column. 
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Figure 9. The stellar mass image of the 
same four galaxies presented in Fig. 8 after 
1.5 Gyrs of simulation, but now the gas is 
replaced by stars. Without gas the bar does, 
or does not, appear just as well. It is obvious, 
though, that gas is needed to generate spiral 
structure. 
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been re-run without gas. In order to obtain nearly the same 
density structure and same rotation curves the gas has been 
replaced by stellar particles which were given the same sech 2 
vertical distribution as the other stars. The resulting density 
structure after 1.5 10 9 years is given in Fig. 9. For v™^c = 
165 and 180 km s _1 the structure closely resembles the mass 
structure in case gas in included. Only the central regions 
are less pronounced because no (gas) mass has been trans- 
ported inwards. Anyway, the conclusion is straightforward; 
gas is not essential for the formation of bars, at least for the 
galactic constitution investigated presently. Yet, when in- 
specting Fig. 9 for the cases of v™^* = 126 and 150 km s _1 
another fact is obvious. Namely that gas is essential for the 
formation of spiral arms. If there is no gas the galactic disc 
is featureless, while if gas is included, the top two rows of 
Fig. 8 show that a rich spiral structure may develop. 



6.2 Q 

Although in principle not applicable for global instabilities, 
Toomre's (1964) Q criterion is often invoked to assess the 
dynamical state of a galactic disc. When rigorously applied, 
Q can then serve as a kind of galactic thermometer, but 
this rigorous application has its problems as we shall see. 
Toomre's Q value indicates when an infinitely thin pure stel- 
lar disc becomes unstable to local instabilities. It is defined 



Q 



3.36Gct ' 



(20) 



where a is the stellar surface density and ft the epicyclic 
frequency given by 



2 _ 2v c dv c 2v c 



(21) 



where v c is the circular test particle rotation. Q is thus de- 
fined there where the epicyclic approximation holds, mean- 
ing (iifj) 1 / 2 < v c - This leads to the first reason why the Q 
criterion does not apply for bars. Bars develop in the central 
regions of galaxies which have a large disc mass contribu- 
tion. At exactly these regions and for those circumstances 
the epicyclic approximation clearly does not hold. The sec- 
ond reason is, of course, that it is a local and not a global 
criterion. 

To give the Q thermometer for the four simulations of 
Fig. 8 it has been calculated as a function of radius. This 
calculation involves a third problem; what to use for v c and 
for a in case of a disc with both stars and gas? If dispersions 
are small the stellar rotation (v t ) is nearly equal to the test- 
particle (~ gas) rotation and it does not matter whether one 
inserts v* or w gas for v c . But in the inner regions of massive 
discs the asymmetric drift is large and it does matter for the 
value of Q which rotation is used. Inclusion of a gas com- 
ponent will decrease the effective Q value as demonstrated 
for a thin disc by Bertin & Romeo (1988). For a real galaxy 
a straightforward recipe for the conversion of Q so as to in- 
clude the gas cannot be given. To make a first attempt one 
might substitute cr* + a gas for the surface density a in Eq. 
(20). In this way I was lead to define two Q values, a pure 
stellar one (Q*) defined as 



Q, 



is 1 / 2 . 



3.36 Gct» ' 
and a Qt given by 



(22) 
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Qt = 



K) 

3.36 G(ct, + a gas ) ' 



(23) 



the two of which should give an unambiguous description 
of the value and range of the Q thermometer. The result 
is given in Fig. 10. The replacement of v* by w gas has an 
effect on Q only at the inner regions. On the other hand, 
for large radii and especially for the least massive discs the 
inclusion of the gas surface density leads to a substantial 
decrease of Q. Only the two most massive discs form a bar. 
Reading the Q thermometer in Fig. 10 it can be noticed 
that for these galaxies both (J* and Qt fall below 1.2 at 
certain inner radii. It is also at these positions where the 
bar eventually develops. Further simulations with different 
galaxy constitutions (next subsection) show that a level of 
1.2± ~ 0.1 for any minimum Q is indeed critical for bar 
(in) stability. 

There is a small complication, however. As noted in 
subsection 4.6, an initially constant vertical to radial velocity 
dispersion ratio of 0.6 changes into a different functionality 
expressed in the function g(R) of Eq. (12) and displayed in 
Fig. 4. For discs with v™^* = 150 km s _1 the setup has been 
adapted to accommodate this changing dispersion ratio. For 
v ™iic = 126 km s _1 this was not necessary because that 
situation settles gently, without mass transport into a stable 
state with a specific g(R). But a different dispersion ratio 
results in a different Q value. If the disc mass is fixed, Q can 
be expressed as 



Q oc 



g(RY 



(24) 



which means that if g(R) is larger than 1.0 as it is in the in- 
ner regions, Q is smaller than when the dispersion ratio were 
0.6 everywhere. As a consequence, for v™* = 126 km s _1 Q 
evolves to slightly lower values a demonstrated in Fig. 10 by 
also showing Q(R) after 1.5 Gyrs. For vjg* = 150 km s _1 
the radial functionality of Q* after 1.5 Gyrs is indistinguish- 
able from the t = situation demonstrating that the input 
g(R) was close to its eventual value. Equation (24) shows 
another important relation. To increase Q by making the 
disc thicker is not easy. For example, to bring up Q from 1.0 
to a stable Q — 2, the value of zo has to be increased by a 
factor four. This changes h/zo = 5 into h/zo = 1.25 and the 
disc is in principle not a disc any more. 
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Figure 10. The initial (t = 0) value of the Q parameters (Q*: 
Eq. 22 and Qf. Eq. 23) for the four galaxies of Fig. 8 (20 = 900 
pc, egN = 25%). For the stable discs Q is also given at t = 1.5 
Gyrs by the dashed lines. When any Q falls below ~ 1.2 at certain 
radii a bar develops. 



6.3 A bar stability criterion 

The matter of bar (in) stability has been further investigated 
by performing a number of simulations for discs with var- 
ious masses and thicknesses. As noted above, the stability 
depends on the precise value of g(R) and as a consequence 
an initial situation should be close to having its equilibrium 
vertical to radial dispersion ratio. For the four different disc 
masses with belonging maximum rotation a good approxi- 
mation to this equilibrium was found by the functionality 
given in Fig. 4. In each case the thickness parameter zo was 
changed until stability changed into instability or vice versa. 

The results of this investigation arc summarized in Ta- 
ble 2 and graphically in the (wJS^, zo) plane in Fig. 11. Even 
for v™c = 126 km s" 1 a bar can be generated when the disc 
is made thin and cool enough. On the other hand, for the 



maximum disc case of w^sc ~ 180 km s _1 it seems impossi- 
ble to avoid a bar. For a very thick disc with zo = 1300 pc a 
bar still appears. To make the disc even thicker, one needs 
large stellar velocity dispersions and the epicyclic approxi- 
mation cannot be applied over most of the radial extent. A 
determination of the proper asymmetric drift then becomes 
problematic. For the more massive discs a long bar develops, 
while the size of the bar decreases towards the less massive 
discs. The dependence of this bar length on the thickness 
of the disc seems to be small. In this way one might deter- 
mine the disc contribution to the total rotation for a barred 
galaxy by relating the length of the bar to the scalelength 
of the disc. Unfortunately, the radial photometric profile of 
barred galaxies is often of an extreme Freeman type II (Free- 
man 1970) making it difficult to assign a scalelength. The 
simulations of barred galaxies also show such an evolution 
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Table 2. Bar stability investigation 
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which has been converted to (v^^,zo) and is plotted as 
the dashed line in Fig. 11. Extrapolation to uJJ** = 180 
km s _1 indicates that a stable disc for that case should be 
achieved approximately near h/zo = 3.1 or zq = 1450 pc, 
which represents a puffed up disc rarely seen in reality (van 
der Kruit & Searle 1982; Kregel et al. 2002). 

Obviously, this criterion is derived for one specific galac- 
tic constitution and its general, dimensionless applicability 
should be further assessed. Does it hold, for example, for 
a different v™*? An indication that it indeed holds comes 
from the study of Bottema & Gerritsen (1997) where a simi- 
lar stability change over has been established for a disc with 
vX* = 120 km s" 1 . The applicability of Eq. (25) will di- 
minish when the galaxy is appreciably different from the 
late type spiral studied in this paper. For example if the 
disc deviates from exponential, when there is a substantial 
bulge, or when the gas content is much larger. It is not ex- 
pected that medium differences in gas content will play an 
important role because as we already saw, the bar develops 
just as well with or without gas. A more general discussion, 
mainly relating to other stability criteria will be given in 
Sect. 9. 

As described above, before the bar develops the disc ex- 
hibits some strong spiral features. These features are tran- 
sient and as a consequence the disc heats up, possibly aided 
by additional scattering on gas accumulations or clouds. 
When a bar is present the surrounding disc is heated by 
emanating arms of variable strength. This heating leads to 
an increased disc thickness which is indicated by the arrows 
towards larger zo values in Fig. 11. As can be seen, when a 
disc is more unstable the heating and thickening is larger, a 
result which does not come as a surprise. Noteworthy is the 
fact that in general when the disc becomes thicker it remains 
approximately equally thick at all radii. This even holds for 
the barred region where the thickness is rigorously calcu- 
lated as the azimuthally averaged rms z-height and using 
zo — 1.10 Zrms- Only when a bulge develops from gas trans- 
ported inwards the thickness remains small or even becomes 
smaller in the bulge region. 



Figure 11. Results of the bar stability investigation for discs 
with various masses and thicknesses. A "U" means that a bar 
develops, in case of an "S" a bar does not develop. The arrows 
indicate the heating and following thickening of the stellar disc 
after 1.5 Gyrs. A division of the (i/jl J* , 20 ) plane in a stable and 
unstable section can be made according to Eq. (25) which is rep- 
resented by the dashed line. 



towards type II and suggest that the original scalelength can 
best be retrieved by taking a global average of the eventual 
profile in the logarithmic domain. 

In Fig. 11a number of sampling points are given along- 
side the instability threshold. Of course, more intermedi- 
ate cases could have been simulated but the results now 
gathered already give a good impression. When the unsta- 
ble (U) and stable (S) data points of Fig. 11 are plotted in 
a \og(zo/h) versus v™J<? hdSt plane, a division by a straight 
line into a stable and unstable section can be made. This 
dimensionless division line is parameterized by 



6.4 Dead or alive 

The present simulations are performed with a rigid, "dead", 
bulge and dark halo potential. What might change when in- 
stead, the bulge and halo would be composed of a collection 
of particles? To answer this question one has to distinguish 
between two aspects: 

A) Bar formation and 

B) Bar persistence 

Let's first consider item A. In first instance the develop- 
ment of the bar depends on the potential of the halo and 
bulge. A more massive halo or bulge increases the potential 
energy and stabilizes against bar formation. For the bulge 
there is an additional effect. When the bulge becomes more 
massive an ILR is forced by the higher circular speed near 
the centre which frustrates the bar formation (Toomre 1981; 
Sellwood & Evans 2001). These effects are the result of a 
deeper, azimuthally symmetric, potential well and conse- 
quently changing from a dead to a life mass distribution 
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will not alter the formation of the bar and the stability cri- 
terion of Eq. (25) remains valid. In second instance, how- 
ever, during the process of bar formation there might be a 
difference if the halo and/or bulge were alive. A life com- 
ponent allows angular momentum transfer more easily and 
resonances between disc and halo might occur. Because the 
maximum disc case (uSl** = 180 km s _1 ) has a relatively 
light halo it is feasible to perform a numerical calculation 
with a life dark halo composed of not too massive particles. 
This has actually been done, using 70k halo plus bulge, 70k 
disc, and 16k gas particles. The result is nearly identical 
to the case with a rigid halo and bulge as presented in the 
bottom row of Fig. 8, both during and at the end of the 
simulation after 1.5 Gyrs. This proves that for at least this 
situation a dynamical interaction between an emerging bar 
with a surrounding life component is not important. 

Bar persistence (item B) is a different story. When an 
already barred galaxy is embedded in a life, non rotating 
isotropic dark halo, the bar is slowed down by dynamical 
friction (Debattista & Sellwood 1998). Eventually the bar 
will then dissolve. This slowing down is severe when the 
galaxy is considerably sub maximum and is nearly absent 
when close to a maximum disc situation. The initial situa- 
tion of Debattista & Sellwood is, of course, rather specific 
and in reality a dark halo will acquire some rotation during 
galaxy formation (Tremaine & Ostriker 1999). According as 
a halo rotates faster, the dynamical friction process is less 
severe and so in reality the bar dissolution will be less strong 
than Debattista & Sellwood claim. Nevertheless, if the sim- 
ulations were to be performed with a life halo there would 
be a dynamical influence on the persistence of the bar. For- 
tunately, as demonstrated above, bars which develope in a 
more sub maximum disc situation tend to be smaller. Such 
bars then reside in a region where the disc density is still 
larger than the dark halo density. Consequently the slow- 
ing down mechanism will only be slightly active if the halo 
were life. The mechanism of Debattista & Sellwood is im- 
portant for sub maximum discs with large bars, a situation 
that maybe never occurs in reality. 

A life bulge may, just as a life halo, also influence the 
persistence of an already existing bar by dynamical fric- 
tion. In addition one might imagine that a small bulge be- 
comes aligned with the bar, enforcing its presence. Concern- 
ing bulges, it should be kept in mind that the investigated 
galaxy has only a small bulge; the surface density of the 
bulge is larger than that of the disc for radii less than 0.9 
kpc. The bars which develope are generally much larger than 
that extent, and any effect of a life bulge on the bar is there- 
fore minor. For an earlier type galaxy with larger bulge the 
situation might be different. 

It would be, nevertheless, interesting to investigate nu- 
merically the effects of a life bulge and halo on the persis- 
tence of the bar. But that is less simple than it might seem. 
At first, realistic and preferably slightly rotating halos and 
bulges have to be created. Secondly, halos are massive and 
large numbers of particles have to be used. If, instead, one 
replaces the halo (and bulge) with less and more massive 
particles inevitably stochastic heating of the bar and sur- 
rounding disc will occur. The bar then dissolves by a nu- 
merical effect and drawing conclusions is difficult. 



7 THE STATE OF THE GAS 
7.1 Rain-clouds 

A good assessment of the state of the gas can be obtained by 
plotting the position for every SPH particle in a temperature 
- density diagram. Two of such diagnostic plots are shown in 
Fig. 12 for two exemplary cases: v™{^ = 126 km s _1 with zo 
= 675 pc and = 165 km s" 1 with z = 900 pc. In all the 
simulations a typical "rain-cloud" structure develops. Most 
of the gas resides around 10 4 K which is the equilibrium 
between FUV heating and the practical upper limit set to 
the temperature by the cooling function. In circumstances of 
high density and/or low FUV flux gas may cool and starts 
to rain down from the cloud. The solid line indicates the 
Jeans instability threshold. If a gas particle falls below this 
line and remains there for a period greater than i spa n, half 
its mass is converted into a stellar particle and the other half 
into a hot SN particle. This rain-cloud structure represents 
a rather uneasy and dynamical equilibrium state. If, for ex- 
ample, the FUV flux is removed, nearly all the gas falls to 
low temperatures within two time steps. Only at low densi- 
ties gas may remain warm by the CR heating. Not shown in 
the diagnostic plot is the hot SN gas. This gas has a tem- 
perature around 1.25 10 7 K (for £sn = 25%) and densities 
in the range of 0.002 to 0.2 M pc~ a . If the hot SN phase 
is switched off the gas particles generally move back to the 
10 4 K phase very quickly. When a bar develops, as for v™^* 
= 165 km s _1 (bottom panel of Fig. 12) a typical extension 
on the 10 4 K cloud towards high density appears. In such 
cases inwards transport of gas generates numerous young 
stars of which the heating outweighs the intrinsic cooling. 

Four temperature ranges have been defined: cold for 
T < 1000 K, luke warm for 1000 K < T < 8000 K, warm 
for 8000 K < T < 10 5 K, and hot (SN) for T > 10 5 K. 
For the top panel in Fig. 12 with = 126 km s _1 the 

total amount of gas is 18.7 10 9 M & of which 13.1 10 9 M & 
resides in the stellar disc region, where R < 25 kpc. Of this 
amount, 18.5%, 25.9%, 54.8%, and 0.9% is in the cold, luke 
warm, warm, and hot phase respectively. In case of v™^ = 
165 km s^ 1 , given in the bottom panel of Fig. 12 the total 
amount of gas is 17.9 10 9 M & of which 12.1 10 9 Mq is at 
R < 25 kpc. Division of that last amount into the cold to 
hot phases gives fractions of 15.8%, 20.6%, 62.6%, and 1.1%. 

There is a considerable amount of luke warm gas; in 
most cases even more than the amount of cold gas. Even if 
one were to increase the cold fraction by increasing the f c pa- 
rameter (Eq. 8) then still the same amount of luke warm gas 
remains. This fact is caused by the abovementioned uneasy 
and dynamical equilibrium state of the gas. SPH particles 
are continuously dropping out of the rain-cloud by natural 
cooling and are pushed up there again by heating. In a man- 
ner, this contradicts the proposed notion of a two phase ISM 
(Field et al. 1969) or three phase ISM (McKee & Ostriker 
1977). Although artificially implemented the hot SN state 
of the gas can be considered as the hot phase of this three 
phase ISM model. 



7.2 The ratio of cold to warm H I 

The whole process of cloud collapse, division in smaller frac- 
tions, further collapse, star formation, bubble formation, etc, 
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Temperature versus density of the gas for two typical 
= 126 km s _1 with zo = 675 pc (top panel) and 
"disc = ^5 km s with zq = 900 pc (bottom panel). Each dot 
represents one gas (SPH) particle. The structure is an equilibrium 
between intrinsic cooling and heating by the FUV field. Once 
a gas particle falls below the Jeans instability line (full drawn 
line) and remains there for a certain time span, it will form a 
stellar particle. The extension at 10 4 K towards high density in 
the bottom panel is typical for gas inside a bar. 



is extremely complex and not well understood. But anyway, 
if at a region gas cools and starts to collapse it takes a certain 
amount of time for the first stars to appear. In the present 
simulations this is implemented by requiring a certain time 
span tspan (Eq. 8) the gas has to remain Jeans unstable and 
thus cold before it may form stars. The longer i sp an, the more 
cold gas is present and one would like to fix this parameter 
by comparing simulated warm to cold gas ratios with actu- 
ally observed ratios in spiral galaxies. Moreover, the effect 
of the amount of SN feedback appears to be slightly related 
to the amount of cold gas and so related to i spa n. 

Unfortunately, observations to determine this ratio are 



difficult and the ensuing analysis is ambiguous. In our local 
Galactic neighbourhood it seems that cold and warm H I are 
evenly divided (Kulkarni & Heiles 1987, Payne et al. 1983), 
though there appears to be luke warm gas (Kalberla et al. 
1985) as well, which considerably complicates the analysis. 
For external galaxies the presence of cold gas has been es- 
tablished (Braun 1997, Young & Lo 1996) yet the amount is 
undetermined. Comparing galaxies, Dickey & Brinks (1993) 
demonstrate that the ratio of cold to warm gas may differ 
considerably from system to system yet its absolute value de- 
pends on the unknown temperature (or temperature range) 
of the cold gas. Finally I came to conclude that, as for now, 
the observations allow a fraction of cold Hi gas anywhere 
between 15 and 80%. This is all rather disappointing and 
means that not much more can be done than quoting gas 
ratios for the various simulations. 

There are, nevertheless, two interesting facts predicted 
by the simulations. At first the presence of a sizable amount 
of luke warm gas. Secondly, a large cold H I gas fraction 
leads to excessive clumping and giant SF regions. If the ratio 
of cold to warm Hi does appear to be large in an actual 
galaxy, an additional stabilization mechanism is needed for 
the cold gas. In this respect magnetic fields might be a good 
candidate. 



7.3 Hi , H 2 , and dust 

An important result of the research by Kennicutt (1989) is 
that the Schmidt law holds for the total (H2 + H 1 ) gas sur- 
face density. This means that the H2 in the inner regions of 
galaxies actively participates in generating instabilities and 
the ensuing star formation. As for most late type galaxies, 
also in the inner regions of NGC 628 there is a gradual in- 
crease of the H2/H1 ratio (see Fig. 2). Following the result 
of Kennicutt, the total gas should and has been included in 
the simulations. But there is a problem. H2 is in a dense 
region where according to the Schmidt law there is a lot of 
SF and consequently a high FUV field. That field should 
quickly dissociate the existing H2 . In real galaxies this does 
not happen and thus another process is needed to maintain 
the H2 which is not incorporated in the simulations. 

Above a certain density threshold the atomic gas be- 
comes molecular. This threshold depends on the pressure, 
the radiation field, and metals plus dust. Especially dust 
is needed to shield the molecules from photodestruction 
(Elmegreen 1993). In addition dust allows H2 to form on 
grain surfaces (Hollenbach & Salpeter 1971). Therefore dust 
explains, at least for a large part, why H2 can exist in the 
presence of the high FUV field. But this also means that 
at positions with large H2 column densities there has to be 
more dust than at other positions. This is in agreement with 
the finding of Huizinga & van Albada (1992) that absorption 
in galaxies mainly takes place in the inner regions. 

Remains the problem that in the simulations the H2 is 
simply replaced with additional H I gas. As noted in the be- 
ginning of this subsection, that has to be done to be able 
to follow the Schmidt law. Yet at the central regions of a 
real galaxy heating and cooling might proceed differently 
compared to the H I mechanism. It will likely proceed at 
lower temperature regimes and FUV flux will be attenu- 
ated by the extra dust. Nevertheless one does not expect 
heating and cooling to be substantially different because the 
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Figure 13. Velocity dispersion of the gas in the three spatial 
directions as a function of radius after 1.5 Gyrs of simulation. The 
dots give the observed z- velocity dispersion in NGC 628 which is, 
inside the stellar disc region (R < 25 kpc), in good agreement 
with the simulated values. 



same metal ions responsible for the cooling (C + and ++ ) 
are present. Compared to dynamical time scales SF is fast 
and in a practical numerical sense always takes place in a 
few timesteps. Considering this, and the presence of only a 
moderate amount of H2 (8.5%), it seems reasonable that the 
replacement of H2 by H I in the simulations will not result in 
errors on scales larger than the employed spatial and time 
resolution. 



7.4 The velocity dispersion of the gas 

If there is no SN feedback but only FUV heating the result- 
ing gas velocity dispersion is in the range of 4 to 7 km s _1 
in the stellar disc region. When SNs are included the gas 
can be stirred up to higher dispersions as was noted by Ger- 
ritsen (1997). For a simulation with v^c = 126 km s~\ 
zq = 675 pc using 40k gas particles (see next section) the 
resulting velocity dispersion at t = 1.5 10 9 years is plotted 
in Fig. 13. By the way, for 16k gas particles the result is 
identical. One can see in Fig. 13 that the velocity disper- 
sion decreases when going to larger radii. This is caused by 
less star formation taking place further out in the disc and 
consequently there is less SN action, being less able to stir 
up the gas. What is interesting is that the radial velocity 
dispersion is significantly larger than the dispersion in the 
other directions, similar to the stellar velocity dispersions. 
The actually observed gas z- velocity dispersion in NGC 628 
(Kamphuis 1993, van der Hulst 1996) is also given in Fig. 
13. There is a good agreement between the observations and 
simulations, which suggests that the employed mechanism of 
SN feedback quite accurately simulates the processes taking 
place in reality. 

Actually matters are more complicated than suggested 
in the paragraph above. At first because the hydrodynamics 



is approximated by a finite resolution numerical scheme. Sec- 
ondly because the H I gas in a real galaxy is not a continuous 
medium, but a hierarchical ensemble of clouds. This ensem- 
ble of clouds cannot be considered as a fluid in equilibrium. 
The clouds and thus the H I can be heated by gravitational 
scattering. In that case a velocity dispersion is predicted of 5 
to 7 km s _1 independent of cloud mass (Jog & Ostriker 1988, 
Gammie et al. 1991). Additional heating can, of course, take 
place through kinetic energy input by SNs. It is unclear if, 
and how accurately the numerical scheme is able to describe 
the detailed process of heating for an ensemble of clouds. 
An indication that there is indeed such a shortcoming is the 
simulated low dispersion values at radii beyond 30 kpc: 2 
to 3 km s _1 is simulated while 5 to 7 km s~ is predicted 
and observed. Obviously the simulated gas is not discrete 
enough or too viscous. 



8 SPIRAL STRUCTURE 
8.1 40k gas simulations 

In order to make a more detailed analysis of spiral structure 
it appeared beneficial to increase the number of gas particles 
from 16.000 to 40.000, giving a gas mass per SPH particle of 
4.7 10 5 Mq. Simulations with the increased number of gas 
particles show lots more details and structure. To illustrate 
this, in Fig. 14 the resulting structure is displayed for four 
cases which have been investigated for various reasons. The 
top two rows are for a galaxy with maximum disc rotation of 
126 km s~\ z = 675 pc, esn = 25% and adding 4 M yr _1 
of new gas. Fig. 15a is after a simulation time of 1.5 Gyrs 
and, except for the larger particle number, is equal to the 
situation in Fig. 5c. Comparing the two one can notice that 
especially in the inner regions more gas particles gives a more 
pronounced spiral structure. In the outer regions the spiral 
arms are more continuous, less flocculent, and contain fewer 
star formation spots. Overall, the general appearance of the 
mass and gas distribution and of the spiral structure has 
not changed by using more gas particles. This implies that 
medium and large scale structures which are generated do 
not depend on the employed particle number. A quantitative 
analysis of structures can thus be made and the result should 
be generally valid. 

In Fig. 14b the simulation of Fig. 14a has been contin- 
ued untill 3 Gyrs to investigate disc heating and the effect of 
an increased young population of stars. The lower two rows 
of Fig. 14 are for a galaxy with maximum disc rotation of 
150 km s~\ z = 900 pc (14c) and z = 675 pc (14d) at 
t = 1.5 Gyrs. These two simulations are on opposite sides 
of the bar instability threshold for that maximum disc ro- 
tation. Fig. 14c is, except for the larger particle number, 
equal to Fig. 8b and the same conclusions apply as for the 
comparison of the = 126 km s _1 discs, given above. 



8.2 Stellar heating 

A long simulation of a stable disc provided the possibility 
to investigate the evolution of the stellar velocity dispersion. 
New stars are formed at a rate of approximately 4 Mq yr _1 
with an original velocity dispersion equal to that of the gas. 
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Figure 14. a to d The resulting structure of the galaxy using 40k gas and 150k stellar particles, a For f^J* = 126 km s _1 , zq = 
pc, at t = 1.5 Gyrs. b As a but now at t = 3 Gyrs. c For = 150 km s _1 , z = 900 pc at t = 1.5 Gyrs. d As c but for z = 675 

More gas particles increases the resolution which results in a more detailed spiral structure specially in the inner regions. 
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Over a time span of a few Gyrs these young stars will acquire 
an ever larger velocity dispersion. 

The proven stable disc with vJJU = 126 km s _1 and zo 
= 675 pc was simulated for 3 Gyrs using 150k stellar, and 
40k gas particles. For the whole period gas was added to 
the system at a rate of 4 MQyr -1 to compensate for the gas 
consumption by star formation. At t = 3 Gyrs, the situation 
displayed in Fig. 14b, the then available stars were divided 
in five age groups. The youngest stars have an age < 0.75 
Gyrs, there are three groups with ages between 0.75 to 1.5 
Gyrs, 1.5 to 2.25 Gyrs, and 2.25 to 3 Gyrs, and there is a 
group of old stars present from the onset of the simulation, 
thus having an age > 3 Gyrs. At radii of 5, 10, 15, and 20 
kpc the velocity dispersion in both the radial and vertical 
direction was determined of these age groups. The result is 
given as a function of the square root of the age of each 
group in Fig. 15. In this way the evolution of the dispersion 
can be monitored at each radial position. As can been seen, 
the increase in dispersion is nicely consistent with a square 
root of age functionality. After 3 Gyrs the dispersion of the 
old population has not yet been reached. If one extrapolates 
the data points to larger ages, then on average, after ~ 6.5 
Gyrs the dispersion of the old, equilibrium population of 
stars will be reached. 

An overview of the rates and the mechanisms of disc 
heating is given by Binney & Tremaine (1987, p. 470 - 484) 
and by Lacey (1991). In short, the collisionless stellar en- 
semble of a stellar disc can be heated by three principal 
mechanisms. Firstly, heating by molecular clouds (Spitzer & 
Schwarzschild 1953, Villumsen 1985), which generates heat- 
ing oc t 0/25 . Secondly, heating by transient spiral features 
(Carlberg & Sellwood 1985) which is more efficient and heats 
oc t ' 5 but molecular clouds are needed in addition to gen- 
erate vertical heating. A third mechanism is heating by the 
infall of small satellite galaxies. Stellar heating as observed 
in the solar neighbourhood indicates that the heating, or 
increase in velocity dispersion, proceeds proportional to i 0,5 
(Wielen 1977) so that the more effective spiral arm heating 
is needed to explain this observation. 

The simulation of Fig. 14a and b actually includes, both 
heating by simulated molecular clouds (conglomerates of gas 
particles), and heating by transient spiral arms. It should 
therefore give a reasonable description of heating as occur- 
ring in real galactic discs. As witnessed in Fig. 15 the sim- 
ulated heating is indeed consistent with a proportionality 
oc t "' , as observed in our Galaxy and predicted by spiral 
arm heating. However, as mentioned in Sect. 4.4 the sim- 
ulation is only collisionless for J?<2/i. Therefore, for radii 
larger than ~ 9 kpc the simulations have additional heating 
by particle scattering. This will generate heating time scales 
being shorter than for real galaxies, certainly at the radial 
positions of 15 and 20 kpc. 



8.3 On what does the amount of spiral structure 
depend? 

My original perception concerning the creation of spiral 
structure was that its origin lies in an underlying young and 
cold stellar fraction. To investigate this assumption, for two 
stable discs (v%£ = 126 km s _1 , z = 675 pc and = 
150 km s _1 with z = 900 pc) long simulations using 40k gas 
particles have been performed until 3 10 9 years. During this 
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Figure 15. The time evolution of the stellar velocity dispersion 
in the radial direction (top) and vertical direction (bottom panel) 
derived from the simulation presented in Fig. 14b. Stellar particles 
are divided up into five age bins at four radial positions of 5, 
10, 15, and 20 kpc. For each bin at each radius the dispersion is 
calculated. The old stars, present from the onset of the simulation 
(age > 3 Gyrs) are indicated by the bars. At each age bin a 
larger dispersion is for a smaller radius. It appears that the stellar 
heating is consistent with being proportional to t 0,5 . 



time young stars have been created at a rate of ~ 4 MQyr -1 
at the position and with the kinematics of the gas layer. The 
result was disappointing. As can be seen in Fig. 14a and b, 
the spiral structure might be marginally stronger at 3 Gyrs 
compared to that after 1.5 Gyrs. However, if the whole time 
sequence for this v™^* = 126 km s^ 1 run is inspected then 
there are no gradual changes any more after ~ 1 Gyrs. There 
are fluctuations in the appearance and strength of the spiral 
arms over periods of ~ 0.2 Gyrs which is illustrated as the 
difference between Fig. 14a and Fig. 14b. The same conclu- 
sion applies for the = 150 km s _1 , zo — 900 pc case 
(Fig. 14c). If the addition of gas is stopped after 1.5 Gyrs, 
the SFR decreases to ~ 2 MQyr -1 at t = 3 Gyrs and the 
spiral structure slowly fades away. 

In retrospect this failure of a young population to gen- 
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erate a more pronounced spiral structure can be understood. 
For example, the wJjU = 126 km s~ x disc has a total stellar 
mass of 45.5 10 9 Mq. During the simulations cold stars are 
added at a rate of 4 Mgyr" 1 . After 1 Gyrs the cold stellar 
fraction then amounts to 4 10 9 Mq, and after 3 Gyrs to 
12 10 9 Mq which is approximately 20% of the then avail- 
able total stellar mass. However, the stars which were formed 
first are then not young and cold any more because of stel- 
lar heating. As can be seen in Fig. 15, after 3 Gyrs, the 
stars with that age have increased their velocity dispersion 
to more than half of that of the old population. Considering 
the numbers and heating rates, the mass fraction of a rea- 
sonably cold stellar population will not surpass the amount 
of 20%. Such a fraction might locally destabilize the disc to 
some extent but will never become a dominant mechanism. 

The total mass of a disc with v™* = 150 km s _1 is 
64.5 10 9 M . In that case for a SFR of 4 M yr _1 the mass 
fraction of cold stars is approximately 10%, being a factor 
two smaller than for the = 126 km s^ 1 disc. For a 

more massive disc, with the same gas content it would thus 
be even more difficult to create instabilities by the cold stel- 
lar population process, and vice versa, of course. Some ad- 
ditional experiments have been carried out by adding large 
quantities of gas during the beginning of the simulations. As 
a result the SFR goes up and a larger cold stellar fraction 
is generated. The spiral structure then becomes more pro- 
nounced, but subsides to its usual state once the extra gas 
flow is stopped. In these cases it is not clear whether the 
more pronounced spiral structure is generated by the larger 
cold stellar fraction or by the larger amount of gas which is 
present. 

Since an increased cold stellar fraction barely enhances 
the spiral structure, the influence of other parameters has 
been investigated. It was apparent from the investigation of 
Sect. 5 that spiral structure is mainly caused by shearing 
gas structures and an increased fraction of cold gas might 
enhance such a process. For the situation of Fig. 14c, v™^* 
= 150 km s _1 and zq — 900 pc, the collapse factor f c (Eq. 
8) has been increased by a factor four in combination with 
or without a change of the supernova efficiency. A larger 
collapse factor produces more cold gas, which is colder and 
has densities up to a factor of 10 larger. The gas arms are 
thinner and more threadlike and tend to form gas clumps 
which cause strong star formation regions. For a combina- 
tion of tspan = 2 x ts and a low SN efficiency of £sn = 
10% the image is like the top rows of Fig. 5 and not like 
any actual galaxy. Consequently, the attempt to make more 
massive gas concentrations by increasing the cold gas frac- 
tion and lowering the SN dispersal mechanism is successful, 
but produces the familiar clumps and does not enhance a 
global spiral pattern. Letting esn remain at 25% or mak- 
ing it larger in combination with a larger cold gas fraction 
does generate the more threadlike gas structure, but does 
not create larger or significantly smaller spiral disturbances 
in the underlying stellar density structure. What this inves- 
tigation makes clear is that the generated spiral structure as 
it is in Fig. 14c is only weakly dependent on the exact value 
of the gas parameters. Putting it in other words, a signifi- 
cant change of cold gas fraction or supernova efficiency does 
not lead to large changes in the spiral structure morphology. 

The standard cooling function applied to the situation 
of NGC 628 leads to a slightly large value of the observable 



SFR, when compared with what is actually observed. Such 
a large SFR had the advantage of quickly generating an ap- 
preciable amount of young and cold stars of which the effect 
on the emerging spiral structure could be investigated. Since 
such effects are demonstrated to be minor, the cooling func- 
tion can be changed so as to generate a lower SFR. This has 
been investigated for the present case (fjjl? = 150 km s _1 , 
z = 900 pc) by lowering the SFR to ~ 2§ M yr _1 and 
adding only 2 Mq yr of gas constantly. It appears that 
there are no noticeable differences in the generated spiral 
structure compared to a situation with a larger SFR. From 
this follows that spiral structure at low gas suppletion rates 
can be sustained longer if the cooling of the gas is less effi- 
cient . 

One can conclude that small and medium scale spiral 
structure as displayed in Fig. 14a to c can be generated in 
an isolated spiral disc, for a large range of gas parameters. 
The principle factor determining the strength of the spiral 
structure is the surface density of the gaslayer. However, it 
is not possible to create autonomously a large scale spiral 
density wave. If the disc is thinner and colder and crosses 
the threshold of bar instability the bar can then drive an 
m — 2 spiral structure as demonstrated in Fig. 14d. 



8.4 A quantitative analysis of spiral structure 

One would like to compare the spiral structure produced 
by the simulations with the actual spiral structure in real 
galaxies in some quantitative way. For example, the number 
of spiral arms, their shape, and strength. Astonishingly only 
a small amount of research has been carried out in this re- 
spect and a comparison of simulations with observations is 
only marginally possible. 

The shapes of spiral arms are generally well repre- 
sented by a logarithmic functionality (Danver 1942, Ken- 
nicutt 1981) as 



roe 



tan(ip) 



(26) 



where r is the position on the arms with belonging an- 
gle O and i p is the pitch angle. Rewriting Eq. (26) gives 
ln(r) oc tan(ip) O. When the density structure of a (face-on) 
galaxy is rebinned to a (ln(r),0) plane, logarithmic spirals 
show up as straight lines. Derivation of the pitch angles for 
each arm is straightforward and the number of crossings of 
one feature over the 2-k long O axis gives the m-mode of 
the feature. On various occasions simulations have been dis- 
played in such a plane and the simulated arms are always 
well represented by a logarithmic shape. However, (ln(r), O) 
planes are not easily comparable because there is in fact 
too much and too dispersed information. One would like to 
reduce the information into less and more informative num- 
bers. 

One method to do this has been explored by Considere 
& Athanassoula (1982, 1988). It involves a Fourier transform 
of the image of the galaxy into a superposition of global 
m-armed logarithmic spirals following Kalnajs (1975). A 
two dimensional grid is transformed into Fourier co- 

efficients A(p,m) by 



A(p,m) 



(27) 
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Figure 16. a to d Amplitude of the Fourier coefficients A(p, m) 
of the V-band images in Fig. 14a to d using global logarithmic 
spirals as base function (Eq. 27). For each panel the lines from 
bottom to top are for m = 1 to 5 and are offset by 0.1 with 
respect to each other. The belonging pitch angle is indicated by 
the numbers for a few features. At panel d the signal at p = 
0, m = 2 originates from the bar of Fig. 14d. 

where p is the wavenumber, m the mode and Wij the value at 
gridpoint i,j. The pitch angle i p is given by tan(i p ) = —m/p 
at the wavenumber where \A(p, m)\ has its largest ampli- 
tude. Considere & Athanassoula (1988) present the calcu- 
lated Fourier coefficients in their Fig. 5 for 16 spiral galax- 
ies, among which are a number known to have a well devel- 
oped grand design structure. For most of the galaxies the 
coefficients |j4(p,m)| rarely surpass the level of 0.1. Notable 



exceptions are NGC 5247 (\A(p,m = 2)| ~ 0.45) and NGC 
6946 (\A(p,m = 2)| ~ 0.21), but these two galaxies both 
have exceptionally large pitch angles of 26 and 32° respec- 
tively. Considere & Athanassoula claim that the two-armed 
component is everywhere dominant, but a close inspection 
of their data shows that only for 8 out of 16 galaxies m = 2 
is dominant, while for 6 out of 16 the m = 1 component 
dominates! Surely this indicates that a number of systems 
have undergone some kind of interaction. NGC 628 is in the 
sample and shows a maximum \A(p, m = 2)| of 0.16 with a 
pitch angle of 16°. In addition there is a peak for m = 3 at 
a level of 0.09 with pitch angle of 28° . 

For the V-band images of the simulated galaxies in Fig. 
14 the Fourier coefficients for modes one to five have been 
calculated following Eq. (27) and are displayed in Fig. 16. 
In general the peak values do not exceed the level of 0.1, 
as for most of the real galaxies, in fact. Fig. 16 is only one 
snapshot; if the Fourier coefficients are monitored over a 
time span, peaks come and go for modes 2 to 5. Time scales 
are typically 0.2 Gyrs and levels are comparable to those dis- 
played in Fig. 16a to c. This behaviour reflects the episodical 
nature of the small and medium scale spiral structure dur- 
ing the simulations. For Fig. 16d, the barred galaxy, peaks 
are more constantly present. The simulations are supposed 
to represent NGC 628. As can be seen in Fig. 16, the values 
of the Fourier coefficients do not reach the observed value. 
In addition, simulated pitch angles are somewhat too large. 
Although the simulated V-band images reasonably resem- 
ble images of actual galaxies, the observed image of NGC 
628 (for example on page 99 of the Shapley-Ames catalog, 
Sandage & Tammann 1981) shows arms which are longer 
and less open. This means that the real galaxy NGC 628 
exhibits more a grand design than its simulated counter- 
part. A further discussion of this matter will follow in Sect. 
10. 

As determined by Kennicutt (1981), pitch angles of late 
type galaxies lie in the range between 12 and 30° and have a 
median value of 18°. As such, the pitch angles of the simu- 
lated galaxies are slightly large but certainly not unusual. 
Note that for the galaxies in common to the samples of 
Considere & Athanassoula and Kennicutt the derived pitch 
angles are in good agreement. But this brings me to the 
question: how representative are both these samples for an 
average galaxy? Probably not representative at all. For both 
samples, galaxies are selected at least partly, for having a 
grand design or otherwise well delineated spiral structure. 
Inevitably this produces a severe bias towards high ampli- 
tude spiral signal and probably towards smaller pitch angles. 

An alternative way to analyse spiral structure has been 
employed by a.o. Elmegreen et al. (1989, 1992), Rix & Zarit- 
sky (1995), and Fuchs & Mollenhoff (1999). According to 
this method, first the underlying smooth galactic distribu- 
tion is removed. Then the residual is divided up in a number 
of circular annuli of which the brightness B is Fourier anal- 
ysed with respect to the azimuthal angle (f>: 

B = ^ c m cos(m(0 - (j> m ))- (28) 

Following Fuchs & Mollenhoff the coefficients c m will be dis- 
played as a function of radius to obtain a direct and quan- 
titative insight into the development of the multiplicity of 
the structure. A determination of the pitch angle of mode m 
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Figure 17. Amplitude of the coefficients of the Fourier sum 
along circular annuli as a function of radius for the four V-band 
images of Fig. 14a to d. Contour levels are at 5, 10, 20, 40, and 
80% of the underlying smooth galactic brightness level. Specific 
features can easily be related to features in Fig. 14 and Fig. 16. 



might be done by plotting the phase <j) m versus In (radius), 
but this has not been investigated presently. 

The V-band images of Fig. 14 have been divided up 
into 35 equally wide annuli ranging from a radius of 0.5 to 
25 kpc. After division by the continuum, the Fourier coeffi- 
cients of Eq. (28) were calculated and are displayed in Fig. 
17. One may notice that a good impression can be obtained 
of the radial functionality of the multiplicity of the spiral 
structure. Comparing Fig. 16 and Fig. 17, certain features 



can be related and put at its proper radius. For example, in 
panel a the small m — 2 signal originates from radii between 
11 and 15 kpc, while the m = 4 signal comes mainly from 
positions in the outer regions. In panel b the m = 3 spiral 
with pitch angle of 29° can obviously be recognized in Fig. 
17. On the other hand, the small m — 2 signal in panel c 
of Fig. 17 at radii < 3 kpc cannot readily be seen in Fig. 
16c. Such a feature appears small in relative units, but if 
this disturbance of ~ 15% is expressed in absolute values 
for NGC 628 it amounts to ~ 25 LqPC -2 , which is certainly 
not negligible. Expressed in relative units, Fig. 17 shows that 
the brightness fluctuations become larger at larger radii. In 
absolute units, however, the signal dominates in the inner 
regions. For the barred galaxy in panel d it can be noticed 
that the m — 2 inner bar drives the m — 2 spiral arms fur- 
ther out, except for the small excursion to m = 3 — 4 around 
a radius of 8 kpc. Of course, during the simulation the spiral 
structure changes and consequently also its representation 
as given in Fig. 17. 

Though this is a nice way to display the spiral structure 
of a galaxy, there is a major problem. There are hardly any 
observations which have been interpreted in this manner and 
hence a comparison between simulations and real data is not 
really possible. Rix & Zaritsky (1995) give ci to c?, values for 
.ff'-band observations of 18 face-on spirals. Their main in- 
tent was to investigate the lopsidedness of galaxies and thus 
they focussed on the c\ coefficient. The near infrared data 
do not extend beyond three scalelengths and the values of 
the C3 to C6 coefficients over that radial extent are relatively 
small, yet consistent with the values in Fig. 17. Only Fuchs 
& Mollenhoff (1999) give a detailed Fourier coefficient anal- 
ysis for one galaxy: NGC 1288. Both, the values and radial 
functionality of the coefficients show a striking similarity 
with the simulations of non-barred galaxies. NGC 1288 is a 
regular, massive galaxy with clear spiral structure but has 
no grand design density wave permeating its disc. As such it 
is comparable with the simulations and the spiral structure 
shows a comparable morphology. 



9 SWING AMPLIFICATION 

The spiral structure which develops in the present numeri- 
cal simulations of an isolated and bar stable galaxy appears 
to be generated by swing amplification. This phenomenon 
where small distortions can be severely amplified in a differ- 
entially rotating (shearing) disc was discoved by Goldreich 
& Lynden-Bell (1965) and the theory has been elaborated 
by o.a. Julian & Toomre (1966), Toomre (1981), Athanas- 
soula (1984), and Athanassoula et al. (1987). As noted by 
Toomre (1981) swing amplification is governed by a conspir- 
acy between three processes: shear, shaking and self-gravity. 
Each of these processes can be specified by a dimensionless 
parameter: F, Q, and X which are the shear rate, the fa- 
miliar Toomre's Q parameter and as it is sometimes called 
Toomre's X parameter. In principle, swing amplification has 
been derived and quantified for a one component, infinitely 
thin disc using the epicyclic approximation. That is a sim- 
plification of reality which should be kept in mind when 
comparing results of the theory with findings for a real or 
numerical galaxy. 

To investigate whether swing amplification is indeed the 
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Figure 18. Swing amplification parameters X(m = 4),T, and 
Q as a function of radius for the galaxies of Fig. 14a and Fig. 14c. 
These parameters are defined in two ways (*: Eq. 22+29+31, t: 
Eq. 23+30+32) depending on the way the influence of the gas is 
included. 



cause of the numerical spiral structure, the theory has been 
applied to the galaxies displayed in Fig. 14. The three above- 
mentioned parameters have to be calculated for the galactic 
constitution which leads again to the problem encountered 
in Sect. 6.2 of what to substitute for the rotation and surface 
density. As in Sect. 6.2 where two kinds of Q are defined also 
here two kinds of F and X will be defined; pure stellar (T* 
and X») given by 



r* = 



.X"* 



r dn* 

Th~dR ' 

2-KGma, ' 



(29) 



(30) 



and total or testparticle parameters (Ft and X t ) given by 

Ft d^lp-aa 
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Here Q and k are the orbital and epicyclic frequencies, a is 
the surface density and m is the azimuthal multiplicity of 
the disturbance which is equal to the number of spiral arms. 
For the galaxies of Fig. 14a (v^ = 126 km s _1 , z = 675 



150 km s" 



zo 



pc, t = 1.5 Gyrs) and Fig. 14c 
900 pc, t = 1.5 Gyrs) the swing amplification parameters 
are given as a function of radius in Fig. 18. As can be seen, 
there is only a marginal difference between T* and F t . On the 
other hand, the Q and X parameter may differ substantially 
mainly depending on whether the gas surface density is, or 
is not included. 

Athanassoula (1984) has calculated the maximum am- 



plifications for a number of different values of the three di- 
mensionless parameters. By interpolating her Fig. 26, the 
amplification can thus be found. But this is a maximum 
amplification provided the arrival phase of the distortion is 
optimum. This means that for an actual case the amplifica- 
tion may be less, or zero, or the amplification of a distortion 
which emerges at a later time with more optimum arrival 
phase may approach this maximum limit. When viewed over 
a longer time span the spiral pattern will thus be rather 
variable just as witnessed in the simulations. Inspection of 
Athanassoula's Fig. 26 establishes some generalities, though. 
For a T value between 0.5 and 1.0 swing amplification is 
only effective for X between typically 0.5 and 2.5 and has 
the largest amplifications when X lies between 1.0 and 1.5. 
A large value of Q dampens the amplification; when Q > 2.0 
the amplification is < 3. When Q is low, say Q < 1.2 then 
the amplifications may exceed the value of 30. 

For the present two galaxies of which the parameters 
are displayed in Fig. 18 the maximum amplifications have 
been derived as a function of radius. To that aim the pure 
stellar (*) and total (t) curves have simply been averaged 
using the assumption that the gas layer and its gravitation 
will influence the swing amplification process to some ex- 
tent. If the amplifications would have been derived for both 
curves separately the emerging picture would not be essen- 
tially different. This is because an increasing Q value results 
in a decreasing amplification with a factor barely dependent 
on r and X. A difference between the X curves mainly ap- 
pears at large values of X where the amplification is small 
anyway. The maximum amplifications for modes 2, 3, 4, 6, 
and 12 are displayed as a function of radius in Fig. 19. 

One may notice that the smaller m-modes (m = 2, 3 ) or 
the two and three armed spiral structures allow the largest 
amplification at the smaller radii. Spiral structures with 
more arms will preferentially develop further out. At cer- 
tain radii the amplifications for a number of different modes 
can approximately be equally large. In practice that means 
that the number of arms can, and likely will, change from 
time to time. For the less massive disc (t>3}jj£ = 126 km s~\ 
top panel of Fig. 19) the amplifications reach at most a fac- 
tor of four which is mainly caused by the relatively large 
value of Q. When that simulation is continued to 3 Gyrs 
as displayed in Fig. 14b the swing amplification parameters 
barely change. Consequently also the amplifications appro- 
priate for that later time are equal to those already given 
for t = 1.5 Gyrs. If the disc is more massive (bottom panel 
of Fig. 19), Q is smaller and the allowed amplifications are 
larger certainly for the smaller m-modes. For the same disc 
but with a smaller thickness (zo = 675 pc instead of 900 
pc) as in Fig. 14d the resulting Q parameter (at t = 0) is 
smaller while the X and F parameters remain nearly equal. 
The emerging amplification scheme for that thinner disc is 
then equal to that in the bottom panel of Fig. 19 except that 
all the amplifications have to be multiplied with a factor 2| 
to 3. That results in values between 30 and 40 for m = 2, 
3, and 4 at radii < 10 kpc. But, as we know a bar devel- 
ops in that case which might be related to these allowed 
large amplifications. When one compares the amplifications 
given in Fig. 19 with the actually emerging modes in Fig. 17 
one may note the exceptional correspondence. For example, 
the m = 3 structure in Fig. 17b between 6 and 11 kpc and 
the same allowed mode by swing amplification at exactly 
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Figure 19. The maximum amplification factors given by swing 
amplification theory for the galaxies of Fig. 14a (top) and Fig. 14c 
(bottom). Note that at certain radii the amplifier may be active 
and equally strong for different arm multiplicities. 



the same radii. Even nicer is the correspondence between 
Fig. 17c and Fig. 19 bottom panel. In that case the radial 
functionality of the appearing structure and the predicted 
amplifications are in striking agreement. 

To investigate if the magnitude of the disturbances is 
consistent with swing amplification theory, for the case of 
Fig. 14c the initial gas distortions have been compared with 
the eventual strength of the stellar mass arms. For example 
at R — 3 kpc gas mass fluctuations are ~ 3.5 Mqpc~ 2 and 
the m = 2 stellar mass arms have a strength of ~ 23.6 
Mqpc" 2 implying an amplification of a factor ~ 6.7. The 
allowed maximum amplification of m = 2 at R — 3 kpc is 
~ 11, in perfect agreement with the simulations. Likewise, 
at R — 5 kpc for m = 3 one has an amplification of ~ 4.7 
while a maximum of ~ 17 is allowed, and at R — 13 kpc 
amplifications are ~ 2.2 while ~ 4.5 is allowed. 

One may conclude that the predictions of swing ampli- 
fication theory and the results of the numerical simulations 



are completely consistent. Not only for the radial positions 
where the modes preferentially occur, but also for the actual 
strength of the spiral arms. 



10 DISCUSSION AND CONCLUSIONS 
10.1 On grand design 

For a late type bar stable galaxy in isolation it does not seem 
possible to have a grand design spiral density wave pervad- 
ing its disc. The spiral structure which does develop is on 
a medium and small scale with a multiplicity increasing as 
a function of radius. That finding can be explained qualita- 
tively and quantitatively by swing amplification; a process 
which occurs naturally in a differentially rotating disc. If 
swing amplification can explain the structure for a late type 
galaxy so well, it is likely that it also applies for other, more 
earlier types. The question which then emerges is whether 
for some galaxy swing amplification allows a spiral structure 
with a constant multiplicity over a sizable range of radii. 
Such a situation would, in principle, be equal to a grand 
design. 

To answer the question one could calculate the ampli- 
fication factors for a number of observed galaxies and in- 
vestigate if a certain mode is dominant over a large radial 
extent. Such a calculation has been done by Athanassoula 
et al. (1987) but unfortunately they do not present the ra- 
dial modal functionality. A global insight can be obtained 
by considering equations 29 and 30 and Fig. 18. To make 
a constant radial multiplicity, the F parameter can only be 
slightly variable and the X parameter needs to be nearly 
constant. For NGC 628 (Fig. 18) that is certainly not the 
case. Making a constant X requires that K 2 R/a is constant, 
while in the mean time the rotation curve should be rather 
flat. That results in the condition that oR ~ constant which 
is impossible to achieve for an exponential disc (a oc e~ R ^ h ). 
It can, however, not be excluded that for some rotation curve 
and photometry a constant X over an extent of, say, two 
scalelengths can be achieved. 

For a number of well known grand design galaxies the 
m = 2 mode is the dominant one. How can such a situation 
arise when swing amplification is the driving mechanism? 
For a late type galaxy (like NGC 628) figures 18 and 19 can 
again give some insight. In case of a maximum disc rotation 
of 126 km s _1 (and zo = 675 pc) the amplifications are given 
in the top panel of Fig. 19. For radii between 3 to 10 kpc the 
m = 4 mode is the largest. Inspection of Eq. (30) shows that 
exactly the same X values are achieved when half the multi- 
plicity (m) is used and a double value of the surface density. 
Then Fig. 19 (top) has also the same appearance but the 
m = 4 is replaced with m = 2, while the amplifications are 
much larger because Q is smaller. As a consequence m = 2 
is dominant between 3 and 10 kpc and might drive a grand 
design m — 2 through the rest of the galaxy. The problem is, 
however, that the surface density has to be doubled, which 
makes the disc massive, with a maximum rotation of ~ 180 
km s _1 and then a bar emerges. So, as soon as a late type 
galaxy reaches the state that m = 2 is active over a con- 
siderable radial extent and with enough amplitude, its disc 
is so massive that instead a bar develops. This explains the 
results of the numerical calculations in this paper, namely 
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that only multiple, non grand design spirals appear, or a 
bar when the disc is more massive. Unless, of course, one 
explicitly prevents the development of a bar (Elmegreen & 
Thomasson 1993). Then a grand design m — 2 may exist in 
a hypothetical galaxy. 

For late type galaxies it seems impossible to create a 
grand design and for other disc like galaxies it seems to be 
very difficult. How does nature then make the well known 
beautiful spiral galaxies which adorn so many front pages 
of astronomy books? What can at first be noticed is that 
these beautiful galaxies all belong to a very small and select 
group; there are simply not many of such galaxies. What 
can further be noticed is that these galaxies generally have 
companions. Kormendy & Norman (1979) already remarked 
that: "in general a galaxy shows a global spiral pattern .... 
by having a density wave driving mechanism in the form of 
a bar or a companion" . For example, well known grand de- 
signs with clear signs of interaction are: M81, M51, M101, 
NGC 7753, NGC 5364, and NGC 4305. An interaction with 
accompanying tidal forcing gives, in principle, a natural ex- 
planation for an emanating m = 2 symmetry. In addition 
both theoretical and numerical calculations (Toomre 1981) 
show that even a modest external forcing may produce a 
transient grand design spiral pattern. Such a pattern when 
viewed from the right angle looks very similar to the ob- 
served morphological and kinematical structure of, for ex- 
ample, M81 (Visser 1980). Considering these arguments and 
the present numerical simulations of isolated galaxies, I then 
conclude that nearly all non barred grand design spirals are 
created by an interaction. 

But how about NGC 628? That galaxy is supposed to be 
the template for a typical late type galaxy, yet its observed 
spiral structure is more of a grand design nature than can 
be generated in the simulations. This apparent contradiction 
can be resolved by looking at the neutral hydrogen obser- 
vations of NGC 628 (Kamphuis & Briggs 1992) which show 
that NGC 628 is not that isolated after all. Instead two gas 
clouds with masses of ~ 10 8 Mq are situated on opposite 
sides of the centre at a radius of ~ 12'. In fact the image of 
these clouds in Fig. 3 of Kamphuis & Briggs is similar to the 
gravitational forcing displayed in Fig. 1 of Toomre (1981). 
That modest forcing already generates a huge transient spi- 
ral wave in the underlying galactic disc. Although the exact 
forcing of the gas clouds in NGC 628 is not known, these 
clouds can certainly explain the grand design structure of 
NGC 628 by swing amplification. 



can be inferred from Fig. 3 of Bottema (1993) such that 
(Vz}](= = 60 ± 10 km s~\ For a scalelength of 4.5 kpc, Eq. 
(33) when related to Fig. 11 implies that < 150 km s _1 , 
even though the thickness of the disc is not known. A more 
massive disc would bring the galaxy into the bar unstable 
region, obviously not consistent with its appearance. 

A different approach to determine the contribution of 
the disc has been employed by Athanassoula et al. (1987). 
They use swing amplification theory to investigate the pres- 
ence of the amplifier for a certain mode in case of differ- 
ent disc contributions. For a sample of 48 galaxies with 
reasonably accurate rotation curves this investigation has 
been performed. It appears that, in order to avoid m = 1, 
for the majority of the sample implies that the disc is less 
massive than the so called "maximum disc situation, hav- 
ing a halo without a hollow core". In practice that means 
Vmtc/VobT « 0.9. The threshold for allowing m = 2 is for 
the majority of the sample at an M/L ratio approximately 
half of that to avoid m = 1 (not surprisingly). In terms of 
rotational values one then has no m = 2 amplification any 
more for /v™ h f < x 0.9 = 0.64. That are the facts 
presented by Athanassoula et al. A problem is, however, that 
it is unclear what is meant with "requiring that the amplifier 
is active". Is a minimum amplification required? And over 
what radial extent is it required to be active? For instance 
can Fig. 19 be regarded as presenting an active m = 2 mode? 

A second problem is associated with the interpretation. 
For example the requirement of no m — 1 is in contradic- 
tion with the results of the Fourier analysis of Considere & 
Athanassoula (1988), which shows that for a large fraction 
of galaxies m = 1 is present and even dominant. On the 
other hand it is not known which fraction of real galaxies 
indeed has m = 2 dominant over a certain radial extent. 
Certainly there are a number of galaxies which do not have 
m — 2 present; which are flocculent or amorphous, which 
implies that v^scA'obT < 0.64 for that category. By using 
the swing amplification theory it is assumed that galaxies 
are not interacting, such that all m = 2 activity has an in- 
ternal origin. That is certainly not valid and a number of 
galaxies could be featureless when they would be isolated. 

Considering the facts presented by Athanassoula et 
al. and the discussion just presented, in my opinion, the 
observed range of morphologies of galaxies is certainly 
consistent with the average (non barred) galaxy having 

max / max _ r> a A 
^disc/^obs — U.D4. 



10.2 On the mass contribution of the disc 

When the stellar velocity dispersion of a disc is measured the 
possibilities for a galaxy's internal constitution are severely 
restrained. For a locally isothermal disc one has the relation 



■ 0.88 (v 



2U/2 



(33) 



Bottema (1993). This equation can be rewritten to zo oc 
( v dSc)~ 2 an d can be drawn onto Fig. 11 once the scale- 
length and dispersion ((t>f )j/Jq) are measured. For example, 
lets take NGC 628 being clearly not barred, which restricts 
the galaxy to lie in the stable section of Fig. 11. Unfortu- 
nately the stellar velocity dispersion is not measured, but 



10.3 On bars 

In Eq. (25) a dimensionless bar stability criterion is given 
for a galaxy comparable to NGC 628. The criterion depends 
on two parameters: v™^* /v™* and the aspect ratio zo/h. 
Here I shall elaborate on other parameters that may be of 
influence for a more general galaxy, but shall not go into 
details on bar formation (see e.g. Lynden-Bell 1979; Toomre 
1981; Binney & Tremaine 1987, p 381). 

Earlier type galaxies have more prominent bulges. 
Adding a bulge to a disc plus halo situation has two ef- 
fects. Firstly, the rotation curve is made flatter in the in- 
ner regions, creating and ILR. Such an ILR may frustrate 
or inhibit the formation of a bar. Secondly adding a bulge 
means that the disc is dynamically less important and the 
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bulge may be considered as a kind of substitute halo. The 
combination of both effects results in a galaxy being more 
stable against bar formation if its bulge is more massive. An 
extreme example of this is presented by Sellwood & Evans 
(2001) where they replace in the inner regions nearly the 
complete dark halo by a bulge. The maximum disc rota- 
tional contribution is then 75% and the minimum Q value 
amounts to 1.75. That galaxy is stable, which is not surpris- 
ing when one considers the rotational contributions and Q 
values presented in Sect. 6. 

The bar forms just as well with or without the amount 
of gas typical for an Sc galaxy. This is contrary to the claim 
by Shlosman & Noguchi (1993) that a considerable gas con- 
tent may inhibit the formation of a bar. Their calculations, 
however, do not include star formation and consequently no 
SF feedback. As a result a few "superclouds" develop which 
rapidly spiral inwards by dynamical friction and strongly 
heat up the stellar content in the inner regions. The galaxy 
is then too hot to form a bar. In reality star formation sets 
in and under normal conditions such massive gas clouds will 
never form. 

Different bar stability criteria have been derived be- 
fore. At first, for example, the T/W criterion by Ostriker 
& Peebles (1973). That criterion is derived for MacLaurin 
spheroids and is as such only applicable for that kind of 
theoretical galaxy. Furthermore, the T and W are not eas- 
ily related to observable quantities. What the calculations 
of Ostriker & Peebles did make very clear was that a large 
spherical (dark) matter contribution is an effective way to 
stabilize against bars. 

A second bar stability criterion put forward by a.o. Sell- 
wood & Carlberg (1984) is the requirement of a specified 
minimum Q value. The problem is to actually determine 
that Q value. As demonstrated in Sect. 6, Q is not unique, 
certainly not when gas is present in a disc. Moreover Q has 
a radial functionality which will be different when the ra- 
dial mass profile is different. One then has to consider the 
minimum Q value, or require Q to be below a certain value 
over a certain radial extent. Even then a small difference in 
Q can mark the change over from stable to unstable. What 
can be determined from the investigations in this paper is 
that stability problems arise as soon as Q (Q» or Q t ) at 
some radius falls below the value of ~ 1.2. It may, however, 
also be 1.1 or 1.3 if the constitution of the galaxy is slightly 
different. 

A third criterion is the one by Efstathiou et al. (1982) 
stating that a disc is stable to bar formation if 



y/GM disc /h ~ 

This criterion ca n be rewri tten considering that for an ex- 
ponential disc ^jGM disc /h = 1.61«dfe£ such that Eq. (34) 
translates into v^^/v^ < 0.56 in order to be stable. Con- 
sidering Fig. 11 which contains the criterion of the present 
study, it is clear that Efstathiou et al's criterion is overly 
restrictive. In addition, for very thin discs with h/zo > 12 it 
is not valid, though such thin discs are probably not present 
in reality. 



10.4 Outlook 

It would be desirable to include the effects of dust and met- 
als in the numerical code. Then matters like FUV shield- 
ing, variable cooling functions, and H2 formation can be in- 
vestigated. That might be important in cases like starburst 
galaxies or ULIRGs. 

A grand design in a non-barred realistic galaxy still has 
to be created. Since that does not appear to be possible in 
an isolated galaxy, encounters have to be simulated. Small 
satellites passing a large galaxy can be simulated using a 
dead halo, but for larger satellites dynamical friction be- 
comes important and a life halo must be used. Re-creation 
of the M51 encounter or the M81-M82 interaction would be 
a nice challenge. 

An obvious extension of the present work is to investi- 
gate bar stability for galaxies with more massive bulges. A 
three parameter stability criterion might then be derived. 

It appears difficult to compare numerical spiral struc- 
ture with actual spiral structure in a quantitative way. This 
is mainly caused by the lack of observational data on spiral 
structure. Surely a program to obtain such data is not that 
difficult. Moreover, asymmetries in galaxies can then be in- 
vestigated by determination of the amount of m = 1 and 
from that the presence and frequency of interactions might 
be established. 

10.5 Conclusions 

Finally a compilation of the main conclusions following from 
the numerical simulations and the subsequent analysis and 
considerations. 

(i) A disc in star formation equilibrium is established by 
two opposing actions. On the one hand gas cooling and dis- 
sipation, on the other hand gas heating by the FUV field of 
young stars and mechanical SN feedback. 

(ii) Too little SN feedback results in a disc with giant SF 
regions; too much SN feedback disperses the gas and dilutes 
the spiral structure. 

(iii) Spiral features are generated in the disc by shear- 
ing gas filements which swing amplify the underlying mass 
distribution. 

(iv) Spiral arms exist for approximately half a rotation 
period and disappear by winding and stretching or by inter- 
action with other spiral features. 

(v) A confirmation is given of Ostriker & Peebles (1973) 
result that a late type spiral galaxy is more unstable to bar 
formation when its disc contribution to the total mass is 
larger. 

(vi) A late type spiral galaxy is more unstable to bar for- 
mation when its thickness or equivalently its stellar velocity 
dispersion is lower. 

(vii) A bar stability criterion is determined for late type 
galaxies, depending on the ratio of maximum disc rotation 
to total rotation and on disc thickness. 

(viii) The bar stability of a disc barely depends on the 
amount of gas. 

(ix) At the maximum disc limit a bar stable galaxy can 
only be created for a disc which is unrealistically thick. 

(x) The simulations display the presence of a substantial 
(20 - 25%) amount of luke warm H I gas in the temperature 
range from 1000 to 8000 K. 
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(xi) In the stellar disc the simulated velocity dispersion 
of the gas is in good agreement with the observations. 

(xii) The increase of the stellar velocity dispersion (stellar 
heating) of newly formed stars proceeds proportional to i 0,5 , 
as observed in the solar neighbourhood. After ~ 6.5 Gyrs 
the velocity dispersion level of the old stellar population is 
reached. 

(xiii) The simulations generate small and medium scale 
spiral structure of which the multiplicity increases as a func- 
tion of radius. 

(xiv) It is not possible to create a grand design m = 2 
spiral. Unless a bar develops which may drive the m = 2 
structure. 

(xv) When compared with some exemplary spiral galax- 
ies, the numerically generated spiral structure is less strong 
and has a larger multiplicity. 

(xvi) The emerging spiral structure is qualitatively and 
quantitatively consistent with the theory of swing amplifi- 
cation. 

(xvii) Swing amplification explains that for a late type 
galaxy a dominant m = 2 mode over a considerable radial 
extent can only exist in a massive disc. Such a disc is un- 
stable to bar formation. As a consequence a grand design 
m — 2 spiral cannot develop in an isolated bar stable late 
type galaxy, just as predicted by the numerical simulations. 

(xviii) The suggestion of Kormendy & Norman (1979) is 
confirmed that most of the grand design m = 2 spirals must 
be generated by an interaction or are driven by a central 
bar. 

(xix) The observed range of morphologies of galaxies is 
consistent with the average (non barred) galaxy having 

..max / max r\ a A 

«di S c/«obs =0.64. 
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